Install packages and load libraries

Install required packages and open their respective libraries: sf, rgdal, and readxl

library(sf)
## Warning: package 'sf' was built under R version 4.2.3
## Linking to GEOS 3.9.3, GDAL 3.5.2, PROJ 8.2.1; sf_use_s2() is TRUE
library(rgdal)
## Warning: package 'rgdal' was built under R version 4.2.3
## Loading required package: sp
## Warning: package 'sp' was built under R version 4.2.3
## The legacy packages maptools, rgdal, and rgeos, underpinning the sp package,
## which was just loaded, will retire in October 2023.
## Please refer to R-spatial evolution reports for details, especially
## https://r-spatial.org/r/2023/05/15/evolution4.html.
## It may be desirable to make the sf package available;
## package maintainers should consider adding sf to Suggests:.
## The sp package is now running under evolution status 2
##      (status 2 uses the sf package in place of rgdal)
## Please note that rgdal will be retired during October 2023,
## plan transition to sf/stars/terra functions using GDAL and PROJ
## at your earliest convenience.
## See https://r-spatial.org/r/2023/05/15/evolution4.html and https://github.com/r-spatial/evolution
## rgdal: version: 1.6-7, (SVN revision 1203)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.5.2, released 2022/09/02
## Path to GDAL shared files: C:/Users/xyaze/AppData/Local/R/win-library/4.2/rgdal/gdal
## GDAL binary built with GEOS: TRUE 
## Loaded PROJ runtime: Rel. 8.2.1, January 1st, 2022, [PJ_VERSION: 821]
## Path to PROJ shared files: C:/Users/xyaze/AppData/Local/R/win-library/4.2/rgdal/proj
## PROJ CDN enabled: FALSE
## Linking to sp version:2.0-0
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading sp or rgdal.
library(readxl)
## Warning: package 'readxl' was built under R version 4.2.3
library(psych)
## Warning: package 'psych' was built under R version 4.2.3
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.2.3
## 
## Attaching package: 'ggplot2'
## The following objects are masked from 'package:psych':
## 
##     %+%, alpha
library(Hmisc)
## Warning: package 'Hmisc' was built under R version 4.2.3
## 
## Attaching package: 'Hmisc'
## The following object is masked from 'package:psych':
## 
##     describe
## The following objects are masked from 'package:base':
## 
##     format.pval, units
library(corrplot)
## Warning: package 'corrplot' was built under R version 4.2.3
## corrplot 0.92 loaded
library(spdep)
## Warning: package 'spdep' was built under R version 4.2.3
## Loading required package: spData
## Warning: package 'spData' was built under R version 4.2.3
library(GWmodel)
## Warning: package 'GWmodel' was built under R version 4.2.3
## Loading required package: robustbase
## Warning: package 'robustbase' was built under R version 4.2.3
## Loading required package: Rcpp
## Warning: package 'Rcpp' was built under R version 4.2.3
## Welcome to GWmodel version 2.3-1.
library(tmap)
## Warning: package 'tmap' was built under R version 4.2.3
## Breaking News: tmap 3.x is retiring. Please test v4, e.g. with
## remotes::install_github('r-tmap/tmap')
library(viridis)
## Warning: package 'viridis' was built under R version 4.2.3
## Loading required package: viridisLite
## Warning: package 'viridisLite' was built under R version 4.2.3
library(RColorBrewer)
library(dplyr)
## Warning: package 'dplyr' was built under R version 4.2.3
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:Hmisc':
## 
##     src, summarize
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(car)
## Loading required package: carData
## Warning: package 'carData' was built under R version 4.2.3
## 
## Attaching package: 'car'
## The following object is masked from 'package:dplyr':
## 
##     recode
## The following object is masked from 'package:psych':
## 
##     logit

Read Nigeria shape file

Read shape file using readOGR() and rename to ngamap NB: Add geopolitical zone map

ngamap1 <- readOGR(dsn = "C:/Users/xyaze/Documents/NGA_adm", layer = "NGA_adm1", verbose = TRUE)
## Warning: OGR support is provided by the sf and terra packages among others
## Warning: OGR support is provided by the sf and terra packages among others
## Warning: OGR support is provided by the sf and terra packages among others
## Warning: OGR support is provided by the sf and terra packages among others
## Warning: OGR support is provided by the sf and terra packages among others
## Warning: OGR support is provided by the sf and terra packages among others
## OGR data source with driver: ESRI Shapefile 
## Source: "C:\Users\xyaze\Documents\NGA_adm", layer: "NGA_adm1"
## with 38 features
## It has 9 fields
## Integer64 fields read as strings:  ID_0 ID_1
class(ngamap1)
## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"
head(ngamap1@data) #view columns
##   ID_0 ISO  NAME_0 ID_1    NAME_1 TYPE_1 ENGTYPE_1 NL_NAME_1 VARNAME_1
## 0  163 NGA Nigeria    1      Abia  State     State      <NA>      <NA>
## 1  163 NGA Nigeria    2   Adamawa  State     State      <NA>      <NA>
## 2  163 NGA Nigeria    3 Akwa Ibom  State     State      <NA>      <NA>
## 3  163 NGA Nigeria    4   Anambra  State     State      <NA>      <NA>
## 4  163 NGA Nigeria    5    Bauchi  State     State      <NA>      <NA>
## 5  163 NGA Nigeria    6   Bayelsa  State     State      <NA>      <NA>
plot(ngamap1) #visualise the Map of Nigeria

geozone <- read_excel("C:/Users/xyaze/Desktop/Geo_Zone.xlsx", sheet = "Sheet1")
ngageo <- sp::merge(ngamap1, geozone, by = 'NAME_1') #merge two data sets
head(ngageo)
## class       : SpatialPolygonsDataFrame 
## features    : 6 
## extent      : 5.376527, 13.72793, 4.270418, 12.5025  (xmin, xmax, ymin, ymax)
## crs         : +proj=longlat +datum=WGS84 +no_defs 
## variables   : 10
## names       :  NAME_1, ID_0, ISO,  NAME_0, ID_1, TYPE_1, ENGTYPE_1, NL_NAME_1, VARNAME_1,    GEO_ZONE 
## min values  :    Abia,  163, NGA, Nigeria,    1,  State,     State,        NA,        NA,  North-East 
## max values  : Bayelsa,  163, NGA, Nigeria,    6,  State,     State,        NA,        NA, South-South
View(ngageo) #view data set in a window
ngageomap <- st_as_sf(ngageo)
class(ngageomap)
## [1] "sf"         "data.frame"
View(ngageomap)
print(ggplot(ngageomap) + geom_sf(aes(fill = GEO_ZONE, color='')) + scale_fill_viridis_d(option = "viridis") + theme_bw()) #plot map

tm_shape(ngageomap) + tm_fill (col = "GEO_ZONE", palette = "viridis", title = "Geopolitical Zones")

##Extract south south data and visualise Plot the map of South-South Nigeria showing the 123 Local Government Areas using the subset() function

ssmap1 <- subset(ngamap1, ID_1 == 3 | ID_1 == 6 | ID_1 == 9 | ID_1 == 10 | ID_1 == 12 | ID_1 == 33) #filter south-south states
ssmap1 <- st_as_sf(ssmap1)
class(ssmap1) #view object type
## [1] "sf"         "data.frame"
head(ssmap1) #snapshot of data -first 6 rows and columns
## Simple feature collection with 6 features and 9 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 4.978159 ymin: 4.270418 xmax: 9.467367 ymax: 7.573423
## Geodetic CRS:  WGS 84
##    ID_0 ISO  NAME_0 ID_1      NAME_1 TYPE_1 ENGTYPE_1 NL_NAME_1 VARNAME_1
## 2   163 NGA Nigeria    3   Akwa Ibom  State     State      <NA>      <NA>
## 5   163 NGA Nigeria    6     Bayelsa  State     State      <NA>      <NA>
## 8   163 NGA Nigeria    9 Cross River  State     State      <NA>      <NA>
## 9   163 NGA Nigeria   10       Delta  State     State      <NA>      <NA>
## 11  163 NGA Nigeria   12         Edo  State     State      <NA>      <NA>
## 32  163 NGA Nigeria   33      Rivers  State     State      <NA>      <NA>
##                          geometry
## 2  MULTIPOLYGON (((7.610695 4....
## 5  MULTIPOLYGON (((6.440416 4....
## 8  MULTIPOLYGON (((8.981203 6....
## 9  MULTIPOLYGON (((5.391528 5....
## 11 MULTIPOLYGON (((6.029233 7....
## 32 MULTIPOLYGON (((6.970418 4....
ggplot(ssmap1) + geom_sf(aes(fill = NAME_1, color='')) + scale_fill_viridis_d(option = "viridis") + theme_bw() #plot map

ngamap <- readOGR(dsn = "C:/Users/xyaze/Documents/NGA_adm", layer = "NGA_adm2", verbose = TRUE)
## Warning: OGR support is provided by the sf and terra packages among others
## Warning: OGR support is provided by the sf and terra packages among others
## Warning: OGR support is provided by the sf and terra packages among others
## Warning: OGR support is provided by the sf and terra packages among others
## Warning: OGR support is provided by the sf and terra packages among others
## Warning: OGR support is provided by the sf and terra packages among others
## OGR data source with driver: ESRI Shapefile 
## Source: "C:\Users\xyaze\Documents\NGA_adm", layer: "NGA_adm2"
## with 775 features
## It has 11 fields
## Integer64 fields read as strings:  ID_0 ID_1 ID_2
ssmap <- subset(ngamap, ID_1 == 3 | ID_1 == 6 | ID_1 == 9 | ID_1 == 10 | ID_1 == 12 | ID_1 == 33) #filter south-south states
ssmap <- st_as_sf(ssmap)
ggplot(ssmap) + geom_sf(aes(fill = NAME_1, color='')) + scale_fill_viridis_d(option = "viridis") + theme_bw()

tm_shape(ssmap) + tm_fill (col = "NAME_1", palette = "viridis", title = "South-South States")

rschdata <- read_excel("C:/Users/xyaze/Desktop/rschdata.xlsx", sheet = "Sheet1") #import data
head(rschdata) #snapshot of data
## # A tibble: 6 × 26
##   OBJECTID NAME_1     NAME_2 LGA   `AREA (km2)`    POP POP_DEN PUBLIC_HF PHF_COV
##      <dbl> <chr>      <chr>  <chr>        <dbl>  <dbl>   <dbl>     <dbl>   <dbl>
## 1      696 Rivers     Obio/… Obio…        263.  708008   2691.        27    3.81
## 2      695 Rivers     Khana  Khana        524.  361163    689.        32    8.86
## 3       69 Akwa Ibom  Uyo    Uyo          154.  408011   2648.        23    5.64
## 4       41 Akwa Ibom  Eket   Eket         210.  314900   1502.        19    6.03
## 5      176 Cross Riv… Calab… Cala…         48.8 190551   3905.        73   38.3 
## 6      118 Bayelsa    Yeneg… Yena…        707.  776979   1099.        45    5.79
## # ℹ 17 more variables: HIV_PREV <dbl>, EST_PLHIV <dbl>, TX_CURR <dbl>,
## #   NOR_TXCURR <dbl>, ART_HF <dbl>, AHF_COV <dbl>, `MPI INDEX` <dbl>,
## #   DEPRIVATION <dbl>, POP_POV <dbl>, `HEALTH_DEP%` <dbl>, HEALTH_DEP <dbl>,
## #   `EDUC_DEP%` <dbl>, EDUC_DEP <dbl>, `WORK_DEP%` <dbl>, WORK_DEP <dbl>,
## #   `SEC_SHOCK%` <dbl>, SEC_SHOCK <dbl>

Merge spatial and non-spatial data

Merge data sets using similar columns

ssdata <- merge(ssmap, rschdata) #merge two data sets
head(ssdata)
## Simple feature collection with 6 features and 35 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 7.465471 ymin: 4.469585 xmax: 8.142541 ymax: 5.164662
## Geodetic CRS:  WGS 84
##      NAME_1        NAME_2 ID_0 ISO  NAME_0 ID_1 ID_2          TYPE_2
## 1 Akwa Ibom          Abak  163 NGA Nigeria    3   39 Local Authority
## 2 Akwa Ibom Eastern Obolo  163 NGA Nigeria    3   40 Local Authority
## 3 Akwa Ibom          Eket  163 NGA Nigeria    3   41 Local Authority
## 4 Akwa Ibom     Esit Eket  163 NGA Nigeria    3   42 Local Authority
## 5 Akwa Ibom      Essien-U  163 NGA Nigeria    3   43 Local Authority
## 6 Akwa Ibom      EtimEkpo  163 NGA Nigeria    3   44 Local Authority
##         ENGTYPE_2 NL_NAME_2   VARNAME_2 OBJECTID           LGA AREA (km2)
## 1 Local Authority      <NA>        <NA>       39          Abak      167.0
## 2 Local Authority      <NA>        <NA>       40 Eastern Obolo      156.5
## 3 Local Authority      <NA>        <NA>       41          Eket      209.7
## 4 Local Authority      <NA>        <NA>       42   Esit - Eket      148.2
## 5 Local Authority      <NA> Essien Udim       43   Essien Udim      298.1
## 6 Local Authority      <NA>   Etim Ekpo       44     Etim Ekpo      200.4
##      POP   POP_DEN PUBLIC_HF   PHF_COV HIV_PREV EST_PLHIV TX_CURR NOR_TXCURR
## 1 154945  927.8144        16 10.326245      5.5  8521.975    5363   3461.228
## 2 125525  802.0767        12  9.559849      5.5  6903.875    5499   4380.801
## 3 314900 1501.6691        19  6.033661      5.5 17319.500   17231   5471.896
## 4  93225  629.0486        15 16.090105      5.5  5127.375    6646   7128.989
## 5 249990  838.6112        21  8.400336      5.5 13749.450    9543   3817.353
## 6 138392  690.5788        14 10.116192      5.5  7611.560    4469   3229.233
##   ART_HF   AHF_COV MPI INDEX DEPRIVATION  POP_POV HEALTH_DEP% HEALTH_DEP
## 1      4  7.458512     0.288        71.3 110475.8        42.0   65076.90
## 2      0  0.000000     0.308        73.2  91884.3        39.4   49456.85
## 3     17  9.865939     0.308        73.2 230506.8        39.4  124070.60
## 4      0  0.000000     0.308        73.2  68240.7        39.4   36730.65
## 5     12 12.574662     0.288        71.3 178242.9        42.0  104995.80
## 6      0  0.000000     0.288        71.3  98673.5        42.0   58124.64
##   EDUC_DEP% EDUC_DEP WORK_DEP% WORK_DEP SEC_SHOCK% SEC_SHOCK
## 1       7.7 11930.77      16.4 25410.98        8.0 12395.600
## 2       6.4  8033.60      16.6 20837.15       15.2 19079.800
## 3       6.4 20153.60      20.0 62980.00        5.6 17634.400
## 4       6.4  5966.40      17.4 16221.15        5.4  5034.150
## 5       7.7 19249.23      16.3 40748.37        5.1 12749.490
## 6       7.7 10656.18      15.8 21865.94        2.7  3736.584
##                         geometry
## 1 MULTIPOLYGON (((7.72106 4.9...
## 2 MULTIPOLYGON (((7.689583 4....
## 3 MULTIPOLYGON (((7.891184 4....
## 4 MULTIPOLYGON (((8.142083 4....
## 5 MULTIPOLYGON (((7.610435 5....
## 6 MULTIPOLYGON (((7.702902 4....
View(ssdata) #view data set in a window
ssdata <- st_as_sf(ssdata) #convert object to sf object
class(ssdata)
## [1] "sf"         "data.frame"
str(ssdata) #view structure of an object
## Classes 'sf' and 'data.frame':   123 obs. of  36 variables:
##  $ NAME_1     : chr  "Akwa Ibom" "Akwa Ibom" "Akwa Ibom" "Akwa Ibom" ...
##  $ NAME_2     : chr  "Abak" "Eastern Obolo" "Eket" "Esit Eket" ...
##  $ ID_0       : chr  "163" "163" "163" "163" ...
##  $ ISO        : chr  "NGA" "NGA" "NGA" "NGA" ...
##  $ NAME_0     : chr  "Nigeria" "Nigeria" "Nigeria" "Nigeria" ...
##  $ ID_1       : chr  "3" "3" "3" "3" ...
##  $ ID_2       : chr  "39" "40" "41" "42" ...
##  $ TYPE_2     : chr  "Local Authority" "Local Authority" "Local Authority" "Local Authority" ...
##  $ ENGTYPE_2  : chr  "Local Authority" "Local Authority" "Local Authority" "Local Authority" ...
##  $ NL_NAME_2  : chr  NA NA NA NA ...
##  $ VARNAME_2  : chr  NA NA NA NA ...
##  $ OBJECTID   : num  39 40 41 42 43 44 45 46 47 48 ...
##  $ LGA        : chr  "Abak" "Eastern Obolo" "Eket" "Esit - Eket" ...
##  $ AREA (km2) : num  167 156 210 148 298 ...
##  $ POP        : num  154945 125525 314900 93225 249990 ...
##  $ POP_DEN    : num  928 802 1502 629 839 ...
##  $ PUBLIC_HF  : num  16 12 19 15 21 14 23 20 18 21 ...
##  $ PHF_COV    : num  10.33 9.56 6.03 16.09 8.4 ...
##  $ HIV_PREV   : num  5.5 5.5 5.5 5.5 5.5 5.5 5.5 5.5 5.5 5.5 ...
##  $ EST_PLHIV  : num  8522 6904 17320 5127 13749 ...
##  $ TX_CURR    : num  5363 5499 17231 6646 9543 ...
##  $ NOR_TXCURR : num  3461 4381 5472 7129 3817 ...
##  $ ART_HF     : num  4 0 17 0 12 0 5 4 0 0 ...
##  $ AHF_COV    : num  7.46 0 9.87 0 12.57 ...
##  $ MPI INDEX  : num  0.288 0.308 0.308 0.308 0.288 0.288 0.283 0.308 0.283 0.283 ...
##  $ DEPRIVATION: num  71.3 73.2 73.2 73.2 71.3 71.3 69.3 73.2 69.3 69.3 ...
##  $ POP_POV    : num  110476 91884 230507 68241 178243 ...
##  $ HEALTH_DEP%: num  42 39.4 39.4 39.4 42 42 40.5 39.4 40.5 40.5 ...
##  $ HEALTH_DEP : num  65077 49457 124071 36731 104996 ...
##  $ EDUC_DEP%  : num  7.7 6.4 6.4 6.4 7.7 7.7 7.1 6.4 7.1 7.1 ...
##  $ EDUC_DEP   : num  11931 8034 20154 5966 19249 ...
##  $ WORK_DEP%  : num  16.4 16.6 20 17.4 16.3 15.8 16.6 16.3 17.4 17.4 ...
##  $ WORK_DEP   : num  25411 20837 62980 16221 40748 ...
##  $ SEC_SHOCK% : num  8 15.2 5.6 5.4 5.1 2.7 15.2 5.1 5.4 5.4 ...
##  $ SEC_SHOCK  : num  12396 19080 17634 5034 12749 ...
##  $ geometry   :sfc_MULTIPOLYGON of length 123; first list element: List of 1
##   ..$ :List of 1
##   .. ..$ : num [1:111, 1:2] 7.72 7.72 7.73 7.73 7.73 ...
##   ..- attr(*, "class")= chr [1:3] "XY" "MULTIPOLYGON" "sfg"
##  - attr(*, "sf_column")= chr "geometry"
##  - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA NA NA NA NA NA NA NA NA ...
##   ..- attr(*, "names")= chr [1:35] "NAME_1" "NAME_2" "ID_0" "ISO" ...

Data Visualisation of the 123 Local Government Areas in South South Nigeria

Map showing the population of South-South Nigeria

describe(rschdata$POP, fast = TRUE)
## rschdata$POP 
##        n  missing distinct     Info     Mean      Gmd      .05      .10 
##      123        0      123        1   246110   140365    97031   121815 
##      .25      .50      .75      .90      .95 
##   154806   215165   284152   442095   526765 
## 
## lowest :  48200  79201  83066  91670  93225, highest: 564498 574199 620174 708008 776979
options(scipen=100000) #before a plot sets legend numbers)
ggplot(rschdata, aes(POP, fill = factor(NAME_1))) + geom_dotplot(stackgroups = TRUE, method = "histodot", stroke=1) +  scale_y_continuous(NULL, breaks = NULL) + scale_fill_viridis_d(option = "viridis") + labs(x = "Population") + theme_bw() + theme(legend.title = element_blank(), legend.position = c(0.9,0.85))
## Bin width defaults to 1/30 of the range of the data. Pick better value with
## `binwidth`.

ggplot(rschdata, aes(POP, fill = factor(NAME_1))) + geom_histogram() +  scale_fill_viridis_d(option = "viridis") + labs(x = "Population", y = "Frequency") + theme_bw() + theme(legend.title = element_blank(), legend.position = c(0.9,0.85))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

tm_shape(ssdata) + tm_fill (col = "POP", palette = "-viridis", title = "Population")

Map showing the population density of South-South Nigeria see if u can change legend scale to 500

describe(rschdata$POP_DEN, fast = TRUE)
## rschdata$POP_DEN 
##        n  missing distinct     Info     Mean      Gmd      .05      .10 
##      123        0      123        1    799.8    805.5     81.1    124.7 
##      .25      .50      .75      .90      .95 
##    268.8    489.9    969.5   1547.8   2643.9 
## 
## lowest : 35.9677 62.8588 73.8839 74.352  74.4682
## highest: 3046.18 3459.24 3904.73 4324.91 5829.43
ggplot(rschdata, aes(POP_DEN, fill = factor(NAME_1))) + geom_dotplot(stackgroups = TRUE, method = "histodot", stroke=1, binwidth = 150) +  scale_y_continuous(NULL, breaks = NULL) + scale_fill_viridis_d(option = "viridis") + labs(x = "Population Density") + theme_bw() + theme(legend.title = element_blank(), legend.position = c(0.9,0.85))

ggplot(rschdata, aes(POP_DEN, fill = factor(NAME_1))) + geom_histogram() +  scale_fill_viridis_d(option = "viridis") + labs(x = "Population Density", y = "Frequency") + theme_bw() + theme(legend.title = element_blank(), legend.position = c(0.9,0.85))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

tm_shape(ssdata) + tm_fill (col = "POP_DEN", palette = "-viridis", title = "Population Density")

Map showing the distribution of Public Health facilities in South-South Nigeria

describe(rschdata$PUBLIC_HF, fast = TRUE)
## rschdata$PUBLIC_HF 
##        n  missing distinct     Info     Mean      Gmd      .05      .10 
##      123        0       46    0.998    28.41    17.69     12.0     14.0 
##      .25      .50      .75      .90      .95 
##     17.0     22.0     34.5     51.6     64.3 
## 
## lowest :   6   8   9  11  12, highest:  65  67  73  92 120
ggplot(rschdata, aes(PUBLIC_HF, fill = factor(NAME_1))) + geom_dotplot(stackgroups = TRUE, method = "histodot", stroke=1, binwidth = 3.5) +  scale_y_continuous(NULL, breaks = NULL) + scale_fill_viridis_d(option = "viridis") + labs(x = "Public Health Facility") + theme_bw() + theme(legend.title = element_blank(), legend.position = c(0.9,0.85))

ggplot(rschdata, aes(PUBLIC_HF, fill = factor(NAME_1))) + geom_histogram() +  scale_fill_viridis_d(option = "viridis") + labs(x = "Public Health Facility", y = "Frequency") + theme_bw() + theme(legend.title = element_blank(), legend.position = c(0.9,0.85))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

tm_shape(ssdata) + tm_fill (col = "PUBLIC_HF", palette = "-viridis", title = "Public HF")+ tm_layout(legend.position = c(0.45, 0.45))

Map showing the coverage of Public Health facilities per 100,000 persons in South-South Nigeria

describe(rschdata$PHF_COV, fast = TRUE)
## rschdata$PHF_COV 
##        n  missing distinct     Info     Mean      Gmd      .05      .10 
##      123        0      123        1    13.99    10.22    3.833    5.615 
##      .25      .50      .75      .90      .95 
##    7.608   10.356   16.622   28.673   33.709 
## 
## lowest : 2.27314 2.78856 3.48071 3.57345 3.64434
## highest: 38.31   38.3735 43.9378 48.1705 72.6141
ggplot(rschdata, aes(PHF_COV, fill = factor(NAME_1))) + geom_dotplot(stackgroups = TRUE, method = "histodot", stroke=1, binwidth = 2.1) +  scale_y_continuous(NULL, breaks = NULL) + scale_fill_viridis_d(option = "viridis") + labs(x = "Public Health Facility Coverage") + theme_bw() + theme(legend.title = element_blank(), legend.position = c(0.9,0.85))

ggplot(rschdata, aes(PHF_COV, fill = factor(NAME_1))) + geom_histogram() +  scale_fill_viridis_d(option = "viridis") + labs(x = "Public Health Facility Coverage", y = "Frequency") + theme_bw() + theme(legend.title = element_blank(), legend.position = c(0.9,0.85))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

tm_shape(ssdata) + tm_fill (col = "PHF_COV", palette = "-viridis", title = "Public HF Coverage")

Map showing the estimated people living with HIV in South-South Nigeria

describe(rschdata$EST_PLHIV, fast = TRUE)
## rschdata$EST_PLHIV 
##        n  missing distinct     Info     Mean      Gmd      .05      .10 
##      123        0      123        1     7578     5389     2347     2909 
##      .25      .50      .75      .90      .95 
##     3717     6057     9634    14703    18687 
## 
## lowest : 964     1578.25 1833.4  1853.79 2014.36
## highest: 20060.3 21450.9 21819.6 22440.6 26904.3
ggplot(rschdata, aes(EST_PLHIV, fill = factor(NAME_1))) + geom_dotplot(stackgroups = TRUE, method = "histodot", stroke=1, binwidth = 850) +  scale_y_continuous(NULL, breaks = NULL) + scale_fill_viridis_d(option = "viridis") + labs(x = "Estimated PLHIV") + theme_bw() + theme(legend.title = element_blank(), legend.position = c(0.9,0.85))

ggplot(rschdata, aes(EST_PLHIV, fill = factor(NAME_1))) + geom_histogram() +  scale_fill_viridis_d(option = "viridis") + labs(x = "Estimated PLHIV", y = "Frequency") + theme_bw() + theme(legend.title = element_blank(), legend.position = c(0.9,0.85))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

tm_shape(ssdata) + tm_fill (col = "EST_PLHIV", palette = "-viridis", title = "Estimated PLHIV")

Map showing the distribution of people currently on HIV treatment (TX_CURR) in South-South Nigeria

describe(rschdata$TX_CURR, fast = TRUE)
## rschdata$TX_CURR 
##        n  missing distinct     Info     Mean      Gmd      .05      .10 
##      123        0      111    0.999     4970     5822        0       15 
##      .25      .50      .75      .90      .95 
##      849     2990     6436    11847    15320 
## 
## lowest :     0    75    96   132   278, highest: 16045 17231 20236 21205 60049
ggplot(rschdata, aes(TX_CURR, fill = factor(NAME_1))) + geom_dotplot(stackgroups = TRUE, method = "histodot", stroke=1) +  scale_fill_viridis_d(option = "viridis") + theme_bw()
## Bin width defaults to 1/30 of the range of the data. Pick better value with
## `binwidth`.

ggplot(rschdata, aes(TX_CURR, fill = factor(NAME_1))) + geom_histogram() +  scale_fill_viridis_d(option = "viridis") + theme_bw()
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

ggplot(ssdata) + geom_sf(aes(fill = TX_CURR)) + scale_fill_viridis(direction = -1) + theme_bw()

Map showing the distribution of people currently on HIV treatment (NORTX_CURR) in South-South Nigeria

describe(rschdata$NOR_TXCURR, fast = TRUE)
## rschdata$NOR_TXCURR 
##        n  missing distinct     Info     Mean      Gmd      .05      .10 
##      123        0      111    0.999     2185     2505    0.000    6.776 
##      .25      .50      .75      .90      .95 
##  399.068 1444.339 3171.242 5384.625 7117.494 
## 
## lowest : 0       33.8798 82.4284 115.571 125.718
## highest: 8420.32 8481.4  10005.2 10895.6 13298.8
ggplot(rschdata, aes(EST_PLHIV, fill = factor(NAME_1))) + geom_dotplot(stackgroups = TRUE, method = "histodot", stroke=1, binwidth = 450) +  scale_y_continuous(NULL, breaks = NULL) + scale_fill_viridis_d(option = "viridis") + labs(x = "Noramlised TX_CURR") + theme_bw() + theme(legend.title = element_blank(), legend.position = c(0.9,0.85))

ggplot(rschdata, aes(EST_PLHIV, fill = factor(NAME_1))) + geom_histogram() +  scale_fill_viridis_d(option = "viridis") + labs(x = "Noramlised TX_CURR", y = "Frequency") + theme_bw() + theme(legend.title = element_blank(), legend.position = c(0.9,0.85))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

ggplot(ssdata) + geom_sf(aes(fill = NOR_TXCURR)) + scale_fill_viridis(direction = -1) + theme_bw()

Map showing the Public Health facilities offering HIV services in South-South Nigeria

describe(rschdata$ART_HF, fast = TRUE)
## rschdata$ART_HF 
##        n  missing distinct     Info     Mean      Gmd      .05      .10 
##      123        0       36    0.997    12.42    11.24      0.0      2.0 
##      .25      .50      .75      .90      .95 
##      5.0     10.0     17.0     27.0     31.8 
## 
## lowest :  0  1  2  3  4, highest: 35 43 47 51 57
ggplot(rschdata, aes(ART_HF, fill = factor(NAME_1))) + geom_dotplot(stackgroups = TRUE, method = "histodot", stroke=1, binwidth = 2.2) +  scale_y_continuous(NULL, breaks = NULL) + scale_fill_viridis_d(option = "viridis") + labs(x = "ART Health Facility") + theme_bw() + theme(legend.title = element_blank(), legend.position = c(0.9,0.85))

ggplot(rschdata, aes(ART_HF, fill = factor(NAME_1))) + geom_histogram() +  scale_fill_viridis_d(option = "viridis") + labs(x = "ART Health Facility", y = "Frequency") + theme_bw() + theme(legend.title = element_blank(), legend.position = c(0.9,0.85))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

tm_shape(ssdata) + tm_fill (col = "ART_HF", palette = "-viridis", title = "ART HF")

Map showing the distribution of Public Health facilities offering HIV services per 10,000 persons on treatment in South-South Nigeria

describe(rschdata$AHF_COV, fast = TRUE)
## rschdata$AHF_COV 
##        n  missing distinct     Info     Mean      Gmd      .05      .10 
##      123        0      105    0.996    97.82    150.2    0.000    0.000 
##      .25      .50      .75      .90      .95 
##    5.154   15.582   74.008  310.402  453.638 
## 
## lowest : 0       1.48251 2.0024  2.57467 2.66449
## highest: 497.238 511.073 569.106 1066.67 1458.33
ggplot(rschdata, aes(AHF_COV, fill = factor(NAME_1))) + geom_dotplot(stackgroups = TRUE, method = "histodot", stroke=1, binwidth = 23) +  scale_y_continuous(NULL, breaks = NULL) + scale_fill_viridis_d(option = "viridis") + labs(x = "ART Health Facility Coverage") + theme_bw() + theme(legend.title = element_blank(), legend.position = c(0.9,0.85))

ggplot(rschdata, aes(AHF_COV, fill = factor(NAME_1))) + geom_histogram() +  scale_fill_viridis_d(option = "viridis") + labs(x = "ART Health Facility Coverage", y = "Frequency") + theme_bw() + theme(legend.title = element_blank(), legend.position = c(0.9,0.85))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

tm_shape(ssdata) + tm_fill (col = "AHF_COV", palette = "-viridis", title = "ART HF Coverage")

Map showing the distribution of population poverty in South-South Nigeria

describe(rschdata$POP_POV, fast = TRUE)
## rschdata$POP_POV 
##        n  missing distinct     Info     Mean      Gmd      .05      .10 
##      123        0      123        1   152156    99149    52314    56797 
##      .25      .50      .75      .90      .95 
##    91804   124811   184881   264225   363202 
## 
## lowest : 31571   34929.3 36120.7 43609.6 47691.1
## highest: 372700  398536  406214  490204  661986
ggplot(rschdata, aes(POP_POV, fill = factor(NAME_1))) + geom_dotplot(stackgroups = TRUE, method = "histodot", stroke=1) +  scale_y_continuous(NULL, breaks = NULL) + scale_fill_viridis_d(option = "viridis") + labs(x = "Population Poverty") + theme_bw() + theme(legend.title = element_blank(), legend.position = c(0.9,0.85))
## Bin width defaults to 1/30 of the range of the data. Pick better value with
## `binwidth`.

ggplot(rschdata, aes(POP_POV, fill = factor(NAME_1))) + geom_histogram() +  scale_fill_viridis_d(option = "viridis") + labs(x = "Population Poverty", y = "Frequency") + theme_bw() + theme(legend.title = element_blank(), legend.position = c(0.9,0.85))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

tm_shape(ssdata) + tm_fill (col = "POP_POV", palette = "-viridis", title = "Population Poverty")

Map showing the distribution of health deprivation in South-South Nigeria

describe(rschdata$HEALTH_DEP, fast = TRUE)
## rschdata$HEALTH_DEP 
##        n  missing distinct     Info     Mean      Gmd      .05      .10 
##      123        0      123        1    88015    48501    36758    44653 
##      .25      .50      .75      .90      .95 
##    55800    74922   104984   160923   186055 
## 
## lowest : 18316   28990   30892.8 31205.2 33939.8
## highest: 199268  200970  214446  235666  247803
ggplot(rschdata, aes(HEALTH_DEP, fill = factor(NAME_1))) + geom_dotplot(stackgroups = TRUE, method = "histodot", stroke=1) +  scale_y_continuous(NULL, breaks = NULL) + scale_fill_viridis_d(option = "viridis") + labs(x = "Health Deprivation") + theme_bw() + theme(legend.title = element_blank(), legend.position = c(0.9,0.85))
## Bin width defaults to 1/30 of the range of the data. Pick better value with
## `binwidth`.

ggplot(rschdata, aes(HEALTH_DEP, fill = factor(NAME_1))) + geom_histogram() +  scale_fill_viridis_d(option = "viridis") + labs(x = "Health Deprivation", y = "Frequency") + theme_bw() + theme(legend.title = element_blank(), legend.position = c(0.9,0.85))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

tm_shape(ssdata) + tm_fill (col = "HEALTH_DEP", palette = "-viridis", title = "Health Deprivation")

Map showing the distribution of education deprivation in South-South Nigeria

attach(rschdata)
plot(EDUC_DEP,TX_CURR, type="p")

describe(rschdata$EDUC_DEP, fast = TRUE)
## rschdata$EDUC_DEP 
##        n  missing distinct     Info     Mean      Gmd      .05      .10 
##      123        0      123        1    21360    14958     7052     8071 
##      .25      .50      .75      .90      .95 
##    11487    17207    26599    41255    54597 
## 
## lowest : 3663.2  5068.86 5430.38 5455.06 5500.2 
## highest: 58282.2 60251.8 60645.8 67416.4 79832.2
ggplot(rschdata, aes(EDUC_DEP, fill = factor(NAME_1))) + geom_dotplot(stackgroups = TRUE, method = "histodot", stroke=1) +  scale_y_continuous(NULL, breaks = NULL) + scale_fill_viridis_d(option = "viridis") + labs(x = "Education Deprivation") + theme_bw() + theme(legend.title = element_blank(), legend.position = c(0.9,0.85))
## Bin width defaults to 1/30 of the range of the data. Pick better value with
## `binwidth`.

ggplot(rschdata, aes(EDUC_DEP, fill = factor(NAME_1))) + geom_histogram() +  scale_fill_viridis_d(option = "viridis") + labs(x = "Education Deprivation", y = "Frequency") + theme_bw() + theme(legend.title = element_blank(), legend.position = c(0.9,0.85))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

tm_shape(ssdata) + tm_fill (col = "EDUC_DEP", palette = "-viridis", title = "Education Deprivation")

Map showing the distribution of work deprivation in South-South Nigeria

describe(rschdata$WORK_DEP, fast = TRUE)
## rschdata$WORK_DEP 
##        n  missing distinct     Info     Mean      Gmd      .05      .10 
##      123        0      123        1    38962    25693    12651    15765 
##      .25      .50      .75      .90      .95 
##    21732    32962    47730    68826    87777 
## 
## lowest : 8531.4  9794.07 9954.53 12125.3 12404.9
## highest: 93594.4 98222.7 98783.4 141602  155396
ggplot(rschdata, aes(WORK_DEP, fill = factor(NAME_1))) + geom_dotplot(stackgroups = TRUE, method = "histodot", stroke=1) +  scale_y_continuous(NULL, breaks = NULL) + scale_fill_viridis_d(option = "viridis") + labs(x = "Work Deprivation") + theme_bw() + theme(legend.title = element_blank(), legend.position = c(0.9,0.85))
## Bin width defaults to 1/30 of the range of the data. Pick better value with
## `binwidth`.

ggplot(rschdata, aes(WORK_DEP, fill = factor(NAME_1))) + geom_histogram() +  scale_fill_viridis_d(option = "viridis") + labs(x = "Work Deprivation", y = "Frequency") + theme_bw() + theme(legend.title = element_blank(), legend.position = c(0.9,0.85))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

tm_shape(ssdata) + tm_fill (col = "WORK_DEP", palette = "-viridis", title = "Work Deprivation")

Map showing the distribution of security shock in South-South Nigeria

describe(rschdata$SEC_SHOCK, fast = TRUE)
## rschdata$SEC_SHOCK 
##        n  missing distinct     Info     Mean      Gmd      .05      .10 
##      123        0      123        1    17555    10894     5666     7366 
##      .25      .50      .75      .90      .95 
##    10583    15130    22254    31170    40245 
## 
## lowest : 2138.43 3736.58 3826.44 4331.96 4410.4 
## highest: 41934.5 43152.6 43510.8 43539.6 49079.7
var(rschdata$SEC_SHOCK)
## [1] 101242585
ggplot(rschdata, aes(SEC_SHOCK, fill = factor(NAME_1))) + geom_dotplot(stackgroups = TRUE, method = "histodot", stroke=1) +  scale_y_continuous(NULL, breaks = NULL) + scale_fill_viridis_d(option = "viridis") + labs(x = "Security Shock") + theme_bw() + theme(legend.title = element_blank(), legend.position = c(0.9,0.85))
## Bin width defaults to 1/30 of the range of the data. Pick better value with
## `binwidth`.

ggplot(rschdata, aes(SEC_SHOCK, fill = factor(NAME_1))) + geom_histogram() +  scale_fill_viridis_d(option = "viridis") + labs(x = "Security Shock", y = "Frequency") + theme_bw() + theme(legend.title = element_blank(), legend.position = c(0.9,0.85))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

tm_shape(ssdata) + tm_fill (col = "SEC_SHOCK", palette = "viridis", title = "SECURITY_SHOCK") 

Aspatial bivariate analysis using Scatter plot and matrix

Scatter plot

rschvar <- rschdata[c('NOR_TXCURR', 'POP_DEN', 'PHF_COV', 'AHF_COV', 'POP_POV', 'HEALTH_DEP', 'EDUC_DEP', 'WORK_DEP', 'SEC_SHOCK')]
C <- cor(rschvar, method = "pearson", use ="complete.obs")
corrplot.mixed(C, tl.cex = 0.8, tl.col="black")

Cc <- cor(rschvar, method = "spearman", use ="complete.obs")
corrplot.mixed(Cc, tl.cex = 0.8, tl.col="black")

pairs(rschvar)

Scatter Plot Matrix

rcorr(as.matrix(rschvar), type ="pearson")
##            NOR_TXCURR POP_DEN PHF_COV AHF_COV POP_POV HEALTH_DEP EDUC_DEP
## NOR_TXCURR       1.00    0.24    0.09   -0.33   -0.06      -0.08    -0.26
## POP_DEN          0.24    1.00   -0.31   -0.20    0.26       0.40     0.15
## PHF_COV          0.09   -0.31    1.00    0.08   -0.33      -0.44    -0.31
## AHF_COV         -0.33   -0.20    0.08    1.00   -0.05      -0.14    -0.04
## POP_POV         -0.06    0.26   -0.33   -0.05    1.00       0.83     0.35
## HEALTH_DEP      -0.08    0.40   -0.44   -0.14    0.83       1.00     0.72
## EDUC_DEP        -0.26    0.15   -0.31   -0.04    0.35       0.72     1.00
## WORK_DEP         0.00    0.41   -0.42   -0.09    0.77       0.89     0.62
## SEC_SHOCK       -0.27    0.20   -0.32   -0.03    0.70       0.78     0.60
##            WORK_DEP SEC_SHOCK
## NOR_TXCURR     0.00     -0.27
## POP_DEN        0.41      0.20
## PHF_COV       -0.42     -0.32
## AHF_COV       -0.09     -0.03
## POP_POV        0.77      0.70
## HEALTH_DEP     0.89      0.78
## EDUC_DEP       0.62      0.60
## WORK_DEP       1.00      0.70
## SEC_SHOCK      0.70      1.00
## 
## n= 123 
## 
## 
## P
##            NOR_TXCURR POP_DEN PHF_COV AHF_COV POP_POV HEALTH_DEP EDUC_DEP
## NOR_TXCURR            0.0087  0.2981  0.0002  0.5160  0.3743     0.0031  
## POP_DEN    0.0087             0.0004  0.0274  0.0034  0.0000     0.0912  
## PHF_COV    0.2981     0.0004          0.3636  0.0002  0.0000     0.0004  
## AHF_COV    0.0002     0.0274  0.3636          0.5887  0.1249     0.6444  
## POP_POV    0.5160     0.0034  0.0002  0.5887          0.0000     0.0000  
## HEALTH_DEP 0.3743     0.0000  0.0000  0.1249  0.0000             0.0000  
## EDUC_DEP   0.0031     0.0912  0.0004  0.6444  0.0000  0.0000             
## WORK_DEP   0.9903     0.0000  0.0000  0.3097  0.0000  0.0000     0.0000  
## SEC_SHOCK  0.0023     0.0246  0.0003  0.7643  0.0000  0.0000     0.0000  
##            WORK_DEP SEC_SHOCK
## NOR_TXCURR 0.9903   0.0023   
## POP_DEN    0.0000   0.0246   
## PHF_COV    0.0000   0.0003   
## AHF_COV    0.3097   0.7643   
## POP_POV    0.0000   0.0000   
## HEALTH_DEP 0.0000   0.0000   
## EDUC_DEP   0.0000   0.0000   
## WORK_DEP            0.0000   
## SEC_SHOCK  0.0000
rcorr(as.matrix(rschvar), type ="spearman")
##            NOR_TXCURR POP_DEN PHF_COV AHF_COV POP_POV HEALTH_DEP EDUC_DEP
## NOR_TXCURR       1.00    0.44   -0.17   -0.32    0.07       0.00    -0.28
## POP_DEN          0.44    1.00   -0.63   -0.30    0.37       0.39     0.05
## PHF_COV         -0.17   -0.63    1.00    0.20   -0.51      -0.62    -0.33
## AHF_COV         -0.32   -0.30    0.20    1.00    0.03       0.02     0.15
## POP_POV          0.07    0.37   -0.51    0.03    1.00       0.84     0.31
## HEALTH_DEP       0.00    0.39   -0.62    0.02    0.84       1.00     0.69
## EDUC_DEP        -0.28    0.05   -0.33    0.15    0.31       0.69     1.00
## WORK_DEP        -0.01    0.41   -0.61   -0.10    0.71       0.86     0.61
## SEC_SHOCK       -0.26    0.12   -0.42    0.09    0.62       0.76     0.63
##            WORK_DEP SEC_SHOCK
## NOR_TXCURR    -0.01     -0.26
## POP_DEN        0.41      0.12
## PHF_COV       -0.61     -0.42
## AHF_COV       -0.10      0.09
## POP_POV        0.71      0.62
## HEALTH_DEP     0.86      0.76
## EDUC_DEP       0.61      0.63
## WORK_DEP       1.00      0.66
## SEC_SHOCK      0.66      1.00
## 
## n= 123 
## 
## 
## P
##            NOR_TXCURR POP_DEN PHF_COV AHF_COV POP_POV HEALTH_DEP EDUC_DEP
## NOR_TXCURR            0.0000  0.0537  0.0003  0.4599  0.9980     0.0020  
## POP_DEN    0.0000             0.0000  0.0006  0.0000  0.0000     0.5998  
## PHF_COV    0.0537     0.0000          0.0303  0.0000  0.0000     0.0002  
## AHF_COV    0.0003     0.0006  0.0303          0.7512  0.8386     0.1093  
## POP_POV    0.4599     0.0000  0.0000  0.7512          0.0000     0.0005  
## HEALTH_DEP 0.9980     0.0000  0.0000  0.8386  0.0000             0.0000  
## EDUC_DEP   0.0020     0.5998  0.0002  0.1093  0.0005  0.0000             
## WORK_DEP   0.8955     0.0000  0.0000  0.2823  0.0000  0.0000     0.0000  
## SEC_SHOCK  0.0038     0.1715  0.0000  0.3035  0.0000  0.0000     0.0000  
##            WORK_DEP SEC_SHOCK
## NOR_TXCURR 0.8955   0.0038   
## POP_DEN    0.0000   0.1715   
## PHF_COV    0.0000   0.0000   
## AHF_COV    0.2823   0.3035   
## POP_POV    0.0000   0.0000   
## HEALTH_DEP 0.0000   0.0000   
## EDUC_DEP   0.0000   0.0000   
## WORK_DEP            0.0000   
## SEC_SHOCK  0.0000

##Near Neighbour Analysis Centroids

ssdatacast <- st_cast(ssdata, "POLYGON") #change multipolygon to polygon in .shf data set
## Warning in st_cast.sf(ssdata, "POLYGON"): repeating attributes for all
## sub-geometries for which they may not be constant
head(ssdatacast)
## Simple feature collection with 6 features and 35 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 7.555695 ymin: 4.472917 xmax: 7.828423 ymax: 5.099989
## Geodetic CRS:  WGS 84
##        NAME_1        NAME_2 ID_0 ISO  NAME_0 ID_1 ID_2          TYPE_2
## 1   Akwa Ibom          Abak  163 NGA Nigeria    3   39 Local Authority
## 2   Akwa Ibom Eastern Obolo  163 NGA Nigeria    3   40 Local Authority
## 2.1 Akwa Ibom Eastern Obolo  163 NGA Nigeria    3   40 Local Authority
## 2.2 Akwa Ibom Eastern Obolo  163 NGA Nigeria    3   40 Local Authority
## 2.3 Akwa Ibom Eastern Obolo  163 NGA Nigeria    3   40 Local Authority
## 2.4 Akwa Ibom Eastern Obolo  163 NGA Nigeria    3   40 Local Authority
##           ENGTYPE_2 NL_NAME_2 VARNAME_2 OBJECTID           LGA AREA (km2)
## 1   Local Authority      <NA>      <NA>       39          Abak      167.0
## 2   Local Authority      <NA>      <NA>       40 Eastern Obolo      156.5
## 2.1 Local Authority      <NA>      <NA>       40 Eastern Obolo      156.5
## 2.2 Local Authority      <NA>      <NA>       40 Eastern Obolo      156.5
## 2.3 Local Authority      <NA>      <NA>       40 Eastern Obolo      156.5
## 2.4 Local Authority      <NA>      <NA>       40 Eastern Obolo      156.5
##        POP  POP_DEN PUBLIC_HF   PHF_COV HIV_PREV EST_PLHIV TX_CURR NOR_TXCURR
## 1   154945 927.8144        16 10.326245      5.5  8521.975    5363   3461.228
## 2   125525 802.0767        12  9.559849      5.5  6903.875    5499   4380.801
## 2.1 125525 802.0767        12  9.559849      5.5  6903.875    5499   4380.801
## 2.2 125525 802.0767        12  9.559849      5.5  6903.875    5499   4380.801
## 2.3 125525 802.0767        12  9.559849      5.5  6903.875    5499   4380.801
## 2.4 125525 802.0767        12  9.559849      5.5  6903.875    5499   4380.801
##     ART_HF  AHF_COV MPI INDEX DEPRIVATION  POP_POV HEALTH_DEP% HEALTH_DEP
## 1        4 7.458512     0.288        71.3 110475.8        42.0   65076.90
## 2        0 0.000000     0.308        73.2  91884.3        39.4   49456.85
## 2.1      0 0.000000     0.308        73.2  91884.3        39.4   49456.85
## 2.2      0 0.000000     0.308        73.2  91884.3        39.4   49456.85
## 2.3      0 0.000000     0.308        73.2  91884.3        39.4   49456.85
## 2.4      0 0.000000     0.308        73.2  91884.3        39.4   49456.85
##     EDUC_DEP% EDUC_DEP WORK_DEP% WORK_DEP SEC_SHOCK% SEC_SHOCK
## 1         7.7 11930.77      16.4 25410.98        8.0   12395.6
## 2         6.4  8033.60      16.6 20837.15       15.2   19079.8
## 2.1       6.4  8033.60      16.6 20837.15       15.2   19079.8
## 2.2       6.4  8033.60      16.6 20837.15       15.2   19079.8
## 2.3       6.4  8033.60      16.6 20837.15       15.2   19079.8
## 2.4       6.4  8033.60      16.6 20837.15       15.2   19079.8
##                           geometry
## 1   POLYGON ((7.72106 4.965561,...
## 2   POLYGON ((7.689583 4.516528...
## 2.1 POLYGON ((7.619028 4.474583...
## 2.2 POLYGON ((7.645416 4.525392...
## 2.3 POLYGON ((7.67736 4.519306,...
## 2.4 POLYGON ((7.555695 4.504077...
st_centroid(ssdata)
## Warning: st_centroid assumes attributes are constant over geometries
## Simple feature collection with 123 features and 35 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 5.236942 ymin: 4.462525 xmax: 9.28852 ymax: 7.342454
## Geodetic CRS:  WGS 84
## First 10 features:
##       NAME_1          NAME_2 ID_0 ISO  NAME_0 ID_1 ID_2          TYPE_2
## 1  Akwa Ibom            Abak  163 NGA Nigeria    3   39 Local Authority
## 2  Akwa Ibom   Eastern Obolo  163 NGA Nigeria    3   40 Local Authority
## 3  Akwa Ibom            Eket  163 NGA Nigeria    3   41 Local Authority
## 4  Akwa Ibom       Esit Eket  163 NGA Nigeria    3   42 Local Authority
## 5  Akwa Ibom        Essien-U  163 NGA Nigeria    3   43 Local Authority
## 6  Akwa Ibom        EtimEkpo  163 NGA Nigeria    3   44 Local Authority
## 7  Akwa Ibom          Etinan  163 NGA Nigeria    3   45 Local Authority
## 8  Akwa Ibom           Ibeno  163 NGA Nigeria    3   46 Local Authority
## 9  Akwa Ibom Ibesikpo Asutan  163 NGA Nigeria    3   47 Local Authority
## 10 Akwa Ibom     Ibiono Ibom  163 NGA Nigeria    3   48 Local Authority
##          ENGTYPE_2 NL_NAME_2   VARNAME_2 OBJECTID             LGA AREA (km2)
## 1  Local Authority      <NA>        <NA>       39            Abak      167.0
## 2  Local Authority      <NA>        <NA>       40   Eastern Obolo      156.5
## 3  Local Authority      <NA>        <NA>       41            Eket      209.7
## 4  Local Authority      <NA>        <NA>       42     Esit - Eket      148.2
## 5  Local Authority      <NA> Essien Udim       43     Essien Udim      298.1
## 6  Local Authority      <NA>   Etim Ekpo       44       Etim Ekpo      200.4
## 7  Local Authority      <NA>        <NA>       45          Etinan      157.2
## 8  Local Authority      <NA>        <NA>       46           Ibeno      271.2
## 9  Local Authority      <NA>        <NA>       47 Ibesikpo Asutan      175.5
## 10 Local Authority      <NA>        <NA>       48     Ibiono Ibom      336.0
##       POP   POP_DEN PUBLIC_HF   PHF_COV HIV_PREV EST_PLHIV TX_CURR NOR_TXCURR
## 1  154945  927.8144        16 10.326245      5.5  8521.975    5363   3461.228
## 2  125525  802.0767        12  9.559849      5.5  6903.875    5499   4380.801
## 3  314900 1501.6691        19  6.033661      5.5 17319.500   17231   5471.896
## 4   93225  629.0486        15 16.090105      5.5  5127.375    6646   7128.989
## 5  249990  838.6112        21  8.400336      5.5 13749.450    9543   3817.353
## 6  138392  690.5788        14 10.116192      5.5  7611.560    4469   3229.233
## 7  275885 1754.9936        23  8.336807      5.5 15173.675    5376   1948.638
## 8  110123  406.0583        20 18.161510      5.5  6056.765   11018  10005.176
## 9  217014 1236.5470        18  8.294396      5.5 11935.770    6594   3038.514
## 10 279438  831.6607        21  7.515084      5.5 15369.090    8008   2865.752
##    ART_HF   AHF_COV MPI INDEX DEPRIVATION   POP_POV HEALTH_DEP% HEALTH_DEP
## 1       4  7.458512     0.288        71.3 110475.79        42.0   65076.90
## 2       0  0.000000     0.308        73.2  91884.30        39.4   49456.85
## 3      17  9.865939     0.308        73.2 230506.80        39.4  124070.60
## 4       0  0.000000     0.308        73.2  68240.70        39.4   36730.65
## 5      12 12.574662     0.288        71.3 178242.87        42.0  104995.80
## 6       0  0.000000     0.288        71.3  98673.50        42.0   58124.64
## 7       5  9.300595     0.283        69.3 191188.30        40.5  111733.43
## 8       4  3.630423     0.308        73.2  80610.04        39.4   43388.46
## 9       0  0.000000     0.283        69.3 150390.70        40.5   87890.67
## 10      0  0.000000     0.283        69.3 193650.53        40.5  113172.39
##    EDUC_DEP%  EDUC_DEP WORK_DEP% WORK_DEP SEC_SHOCK% SEC_SHOCK
## 1        7.7 11930.765      16.4 25410.98        8.0 12395.600
## 2        6.4  8033.600      16.6 20837.15       15.2 19079.800
## 3        6.4 20153.600      20.0 62980.00        5.6 17634.400
## 4        6.4  5966.400      17.4 16221.15        5.4  5034.150
## 5        7.7 19249.230      16.3 40748.37        5.1 12749.490
## 6        7.7 10656.184      15.8 21865.94        2.7  3736.584
## 7        7.1 19587.835      16.6 45796.91       15.2 41934.520
## 8        6.4  7047.872      16.3 17950.05        5.1  5616.273
## 9        7.1 15407.994      17.4 37760.44        5.4 11718.756
## 10       7.1 19840.098      17.4 48622.21        5.4 15089.652
##                     geometry
## 1  POINT (7.760623 5.006096)
## 2  POINT (7.696764 4.516807)
## 3   POINT (7.94293 4.628248)
## 4  POINT (8.067648 4.633393)
## 5  POINT (7.638981 5.099096)
## 6  POINT (7.590369 4.960887)
## 7  POINT (7.834807 4.801225)
## 8   POINT (8.09915 4.561915)
## 9  POINT (7.940853 4.911513)
## 10 POINT (7.881161 5.199924)
plot(st_geometry(ssdata))
plot(st_geometry(st_centroid(ssdata)), pch = 3, cex = 0.5, col = 'red', add=TRUE)
## Warning: st_centroid assumes attributes are constant over geometries

Near Neighbour Analysis

ID_2 = row.names(ssdatacast)
ssdatacast$ID_2 = ID_2
sscastnb1 <- poly2nb(as(ssdata, "Spatial")) # creating adjacency matrix (neighborhood graph)
summary(sscastnb1)
## Neighbour list object:
## Number of regions: 123 
## Number of nonzero links: 630 
## Percentage nonzero weights: 4.164188 
## Average number of links: 5.121951 
## Link number distribution:
## 
##  1  2  3  4  5  6  7  8  9 10 11 
##  2  5 12 34 24 18 17  5  2  3  1 
## 2 least connected regions:
## 43 119 with 1 link
## 1 most connected region:
## 95 with 11 links
head(sscastnb1)
## [[1]]
## [1]  5  6 12 18 26 28 31
## 
## [[2]]
## [1]   8  13  18  24 120
## 
## [[3]]
## [1]  4  7  8 21 24
## 
## [[4]]
## [1]  3  8 17 21 23 30
## 
## [[5]]
## [1]  1  6 11 12 14 22
## 
## [[6]]
## [1]  1  5 11 28
plot(st_geometry(ssdata), border="grey60", axes = TRUE)
plot(sscastnb1, st_coordinates(st_centroid(ssdata)), pch = 19, cex = 0.01, col="red", add=TRUE)
## Warning: st_centroid assumes attributes are constant over geometries

head(st_coordinates(st_centroid(ssdata)))
## Warning: st_centroid assumes attributes are constant over geometries
##             X        Y
## [1,] 7.760623 5.006096
## [2,] 7.696764 4.516807
## [3,] 7.942930 4.628248
## [4,] 8.067648 4.633393
## [5,] 7.638981 5.099096
## [6,] 7.590369 4.960887

Near Neighbor 2

sscastcoords <- st_coordinates(st_centroid(st_geometry(ssdata)))
sscastnbsf = nb2lines(sscastnb1, coords=sscastcoords, proj4string=CRS("+proj=longlat +datum=WGS84"), as_sf=T) # weights are binary
summary(sscastnbsf)
##        i                j              i_ID               j_ID          
##  Min.   :  1.00   Min.   :  1.00   Length:630         Length:630        
##  1st Qu.: 30.25   1st Qu.: 30.25   Class :character   Class :character  
##  Median : 63.00   Median : 63.00   Mode  :character   Mode  :character  
##  Mean   : 61.89   Mean   : 61.89                                        
##  3rd Qu.: 93.00   3rd Qu.: 93.00                                        
##  Max.   :123.00   Max.   :123.00                                        
##        wt             geometry  
##  Min.   :1   LINESTRING   :630  
##  1st Qu.:1   epsg:NA      :  0  
##  Median :1   +proj=long...:  0  
##  Mean   :1                      
##  3rd Qu.:1                      
##  Max.   :1
ggplot(ssdata) + geom_sf(fill = 'NA', color = 'grey', border='grey60') + geom_sf(data = sscastnbsf, pch = 3, cex = 0.5, col="red") + theme_minimal() + ylab("Latitude") + xlab("Longitude")
## Warning in layer_sf(geom = GeomSf, data = data, mapping = mapping, stat = stat,
## : Ignoring unknown parameters: `border`

Moran I Test

moran(ssdata$NOR_TXCURR, nb2listw(sscastnb1, style="W", zero.policy=TRUE), length(sscastnb1), Szero(nb2listw(sscastnb1, style="W", zero.policy=TRUE)), zero.policy = TRUE)
## $I
## [1] 0.3333579
## 
## $K
## [1] 6.833905
moran.test(ssdata$NOR_TXCURR, nb2listw(sscastnb1, style="W", zero.policy=TRUE), na.action=na.pass, zero.policy = TRUE, alternative = "two.sided") #zero.policy for missing values, na.pass puts a 0 instead of deleting observation)
## 
##  Moran I test under randomisation
## 
## data:  ssdata$NOR_TXCURR  
## weights: nb2listw(sscastnb1, style = "W", zero.policy = TRUE)    
## 
## Moran I statistic standard deviate = 6.024, p-value = 0.000000001701
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.333357854      -0.008196721       0.003214733
moran.mc(ssdata$NOR_TXCURR, nb2listw(sscastnb1, style = "W", zero.policy = TRUE), zero.policy = TRUE, nsim = 999, alternative = "two.sided")
## 
##  Monte-Carlo simulation of Moran I
## 
## data:  ssdata$NOR_TXCURR 
## weights: nb2listw(sscastnb1, style = "W", zero.policy = TRUE)  
## number of simulations + 1: 1000 
## 
## statistic = 0.33336, observed rank = 1000, p-value <
## 0.00000000000000022
## alternative hypothesis: two.sided
moran.plot(ssdata$NOR_TXCURR, nb2listw(sscastnb1, zero.policy=TRUE))

moran.plot(as.vector(scale(ssdata$NOR_TXCURR)), nb2listw(sscastnb1, zero.policy=TRUE), xlim=c(-0.8, 8), ylim=c(-0.8,2.8), pch=1, ylab = 'Spatially Lagged NOR_TXCURR', xlab = 'Normalised TX_CURR')

sslmoran <-localmoran(ssdata$NOR_TXCURR, nb2listw(sscastnb1, style="W", zero.policy=TRUE), zero.policy = TRUE, na.action=na.pass)
print(sslmoran)
##               Ii            E.Ii         Var.Ii       Z.Ii  Pr(z != E(Ii))
## 1    0.262799440 -0.002141948541  0.03569417432  1.4023330 0.1608158453402
## 2    0.453875507 -0.006340871331  0.14987250485  1.1887789 0.2345266842014
## 3    1.606879977 -0.014208229516  0.33316605836  2.8085119 0.0049771040414
## 4    3.551788973 -0.032145961061  0.61145255071  4.5833048 0.0000045768422
## 5    0.339062502 -0.003504176532  0.06862588206  1.3076790 0.1909822112521
## 6    0.153755786 -0.001433974465  0.04293979441  0.7489160 0.4539078643927
## 7   -0.041213007 -0.000073497388  0.00144433011 -1.0824950 0.2790326542547
## 8    4.543410962 -0.080427934843  1.45351382586  3.8352450 0.0001254392000
## 9    0.129015453 -0.000957989160  0.02276564423  0.8614194 0.3890070965404
## 10   0.088040343 -0.000609407414  0.01448700431  0.7365262 0.4614105042275
## 11   0.110555948 -0.000348649238  0.02125730822  0.7606685 0.4468551106494
## 12   0.156166229 -0.000444631917  0.00742210394  1.8178513 0.0690868588962
## 13  -0.366147633 -0.004179024541  0.09899010613 -1.1504692 0.2499506550092
## 14   0.448276411 -0.030668534366  1.19870075876  0.4374518 0.6617837272740
## 15  -0.030969898 -0.001322395264  0.03960302440 -0.1489786 0.8815705047423
## 16   0.289536895 -0.001975592699  0.05912624955  1.1988559 0.2305839685432
## 17   6.280010348 -0.099786205220  2.13674043913  4.3644643 0.0000127434661
## 18  -0.051854263 -0.000028272201  0.00040952644 -2.5609826 0.0104376569128
## 19   0.034514414 -0.000008231447  0.00024683933  2.1973381 0.0279963073197
## 20  -0.283654201 -0.005814270396  0.17334228101 -0.6673327 0.5045596417553
## 21  -0.008932743 -0.000001068861  0.00001548301 -2.2698919 0.0232141412705
## 22  -0.154286451 -0.000300587030  0.00901117517 -1.6221467 0.1047719322864
## 23   1.861093252 -0.048346784180  0.58722155218  2.4917532 0.0127114316274
## 24   0.943028028 -0.009841586946  0.19151213394  2.1773855 0.0294518151042
## 25   2.931836118 -0.027463530972  0.63532781279  3.7127013 0.0002050588576
## 26   0.461585200 -0.010686159535  0.20776970272  1.0360972 0.3001567987920
## 27   3.364479804 -0.022328773206  0.65463535048  4.1859200 0.0000284013372
## 28   0.087935644 -0.000080568057  0.00241584829  1.7907207 0.0733381308427
## 29   1.006944880 -0.020916896615  0.40247905588  1.6201793 0.1051937713384
## 30   0.923285609 -0.001349087713  0.03204715556  5.1650608 0.0000002403607
## 31   0.399350353 -0.010124858943  0.11410565364  1.2121993 0.2254360931278
## 32   0.419899345 -0.002618498178  0.06212253706  1.6951981 0.0900378429644
## 33   0.511401589 -0.005565571029  0.16596925275  1.2689635 0.2044540765636
## 34   0.465561924 -0.006279078720  0.25159720307  0.9406819 0.3468679126656
## 35   0.455194849 -0.003935072872  0.15804686124  1.1548954 0.2481332545047
## 36   0.362016241 -0.003076381111  0.06027379011  1.4870953 0.1369896382763
## 37   0.290721050 -0.004538262499  0.06544082860  1.1541951 0.2484202099147
## 38   0.495395881 -0.005155176712  0.08564823995  1.7103664 0.0871981361078
## 39   0.031060297 -0.000028177181  0.00055374750  1.3211238 0.1864600830176
## 40   0.222690004 -0.004025627438  0.24454159173  0.4584639 0.6466192282599
## 41   0.064132501 -0.000182334839  0.00304445813  1.1656178 0.2437690159198
## 42  -0.856572069 -0.003895904060  0.09230996323 -2.8064663 0.0050088148951
## 43  -3.069109149 -0.162441410889 16.73466646684 -0.7105373 0.4773710307771
## 44   0.131445909 -0.006279078720  0.25159720307  0.2745743 0.7836433366700
## 45   0.248690770 -0.004169508458  0.12451223711  0.7165964 0.4736231956126
## 46   0.117754167 -0.003280919562  0.07778645034  0.4339696 0.6643105164214
## 47  -1.063438040 -0.051131672413  1.95632676267 -0.7237545 0.4692164739028
## 48   0.037642169 -0.000009794150  0.00023296898  2.4668263 0.0136316451239
## 49   0.113495716 -0.006279078720  0.25159720307  0.2387880 0.8112699620346
## 50  -0.156529041 -0.001264667477  0.02482294154 -0.9854741 0.3243912701523
## 51   0.288744767 -0.003164862206  0.19241952395  0.6654629 0.5057544601694
## 52   0.043856871 -0.002758143368  0.08248198284  0.1623104 0.8710614460200
## 53   0.123082893 -0.000721542122  0.02162170668  0.8419595 0.3998106071470
## 54  -0.384187684 -0.003176207340  0.03604665784 -2.0068068 0.0447702424843
## 55  -0.191306225 -0.001809947673  0.04297489807 -0.9140990 0.3606648056928
## 56  -0.041627846 -0.000049474610  0.00148355158 -1.0794838 0.2803720900510
## 57   0.092412711 -0.003189152856  0.09533005594  0.3096360 0.7568377846068
## 58   0.572555629 -0.004641241212  0.13853373374  1.5507657 0.1209578337871
## 59   0.222094682 -0.004180346625  0.06952044265  0.8581839 0.3907909232812
## 60   0.056295509 -0.000040286231  0.00095824098  1.8198980 0.0687745340556
## 61   0.290816374 -0.006279078720  0.18711220564  0.6868232 0.4921941593385
## 62   0.121474067 -0.000594351236  0.01412929747  1.0269352 0.3044509877977
## 63   0.209191113 -0.000832308087  0.02493815123  1.3299509 0.1835344562680
## 64   0.433641344 -0.004303022046  0.08420293704  1.5092296 0.1312401097096
## 65  -0.183487935 -0.000588773016  0.01764549633 -1.3768758 0.1685506406435
## 66   0.300408454 -0.003508015150  0.10482793414  0.9386760 0.3478971116224
## 67   0.296058013 -0.002451223516  0.07332613800  1.1023725 0.2702997645785
## 68   0.074923887 -0.000339590729  0.00491748363  1.0732795 0.2831457464410
## 69   0.271924710 -0.003025616239  0.05037526962  1.2250267 0.2205651680142
## 70   0.334743568 -0.005204003728  0.08645520982  1.1561558 0.2476174385511
## 71   0.325323129 -0.004986583678  0.14879002095  0.8563168 0.3918225963041
## 72  -0.352626141 -0.003182789262  0.12792895718 -0.9769950 0.3285716145660
## 73   0.384928186 -0.005632420217  0.11006987794  1.1772107 0.2391114269497
## 74  -0.197314166 -0.001062137328  0.03181712275 -1.1002306 0.2712316804429
## 75  -0.250683659 -0.002109326935  0.06312023673 -0.9894001 0.3224674024121
## 76   0.311845498 -0.001879438618  0.02135742610  2.1467160 0.0318159016602
## 77   0.223297680 -0.004521422177  0.07516687934  0.8309538 0.4059997323292
## 78   0.345017813 -0.003549084637  0.10605081830  1.0703582 0.2844581187779
## 79   0.318080212 -0.003640835827  0.08628843937  1.0952245 0.2734182764771
## 80   0.105654545 -0.000517322521  0.00863487784  1.1425662 0.2532187625257
## 81   0.078696762 -0.006279078720  0.12262720822  0.2426622 0.8082671268726
## 82  -0.137694467 -0.000776302934  0.01845140965 -1.0079669 0.3134703310049
## 83   0.565098123 -0.005135797196  0.20602366942  1.2563033 0.2090060386896
## 84   0.089015647 -0.000177007017  0.00530707665  1.2243382 0.2208246627724
## 85   0.287257549 -0.001944426205  0.05819530477  1.1988291 0.2305943995008
## 86   0.142726903 -0.000333170293  0.00792240702  1.6072738 0.1079943173271
## 87   0.578669726 -0.005577320481  0.10899914742  1.7696402 0.0767870976051
## 88   0.564855511 -0.005312316044  0.10384775806  1.7693109 0.0768420109934
## 89   0.520546970 -0.006279078720  0.18711220564  1.2179128 0.2232571261932
## 90   0.625521264 -0.006279078720  0.25159720307  1.2595835 0.2078196617412
## 91   0.328960688 -0.001866956745  0.03112018934  1.8753429 0.0607455715388
## 92   0.576489111 -0.006279078720  0.10420292324  1.8053270 0.0710235207650
## 93   0.363532456 -0.006279078720  0.10420292324  1.1456198 0.2519525106859
## 94  -0.128226605 -0.000336903287  0.01358014341 -1.0974466 0.2724462122745
## 95   0.390899548 -0.003928033488  0.04013432229  1.9708316 0.0487431412103
## 96   0.421166472 -0.006279078720  0.10420292324  1.3241612 0.1854495224579
## 97   0.403632792 -0.004209567905  0.25566805064  0.8065924 0.4199013806112
## 98   0.548449695 -0.005105481269  0.15231949193  1.4183492 0.1560888325270
## 99   0.590218568 -0.004083741089  0.12196150629  1.7017511 0.0888020455516
## 100  0.622149899 -0.006279078720  0.10420292324  1.9467772 0.0515614513102
## 101  0.388036777 -0.004664206378  0.07752948910  1.4103553 0.1584348038729
## 102  0.002580034 -0.000000453478  0.00001359871  0.6997660 0.4840734551855
## 103  0.123357333 -0.003805859225  0.09018458304  0.4234433 0.6719718542450
## 104  0.491844751 -0.006085869673  0.24390289803  1.0082317 0.3133432364167
## 105 -0.210218418 -0.006279078720  0.18711220564 -0.4714655 0.6373083277561
## 106  0.358491796 -0.005535466117  0.33574881421  0.6282416 0.5298457064391
## 107 -0.165031090 -0.002566726942  0.05031413386 -0.7242910 0.4688870768594
## 108  0.020238390 -0.001265094420  0.01830238870  0.1589480 0.8737098788773
## 109 -0.270436816 -0.000990209999  0.01944122076 -1.9324620 0.0533025136999
## 110 -0.071160894 -0.001885458427  0.03142801018 -0.3907696 0.6959675118241
## 111 -0.138727146 -0.000397493033  0.01191512533 -1.2672607 0.2050620754513
## 112  0.229389483 -0.006922429959  0.16352243355  0.5843822 0.5589632305307
## 113 -0.021213163 -0.000009997125  0.00040310319 -1.0560697 0.2909363503425
## 114  0.086941147 -0.017871166983  0.22401497393  0.2214490 0.8247428182119
## 115 -0.348641590 -0.052138394127  0.82531953386 -0.3263761 0.7441398363457
## 116 -0.191330838 -0.002284054991  0.04478576003 -0.8933041 0.3716943550869
## 117  0.060382436 -0.000359717976  0.00855345328  0.6567793 0.5113228454889
## 118 -0.084801368 -0.001533515921  0.03642146236 -0.4363135 0.6626092724731
## 119  0.080282714 -0.001089418899  0.13385254408  0.2224141 0.8239915098185
## 120 -0.167793089 -0.006279078720  0.18711220564 -0.3733870 0.7088604661011
## 121  0.351800804 -0.001293655837  0.03073209448  2.0141640 0.0439923206612
## 122  0.061032925 -0.000099646386  0.00298785855  1.1183884 0.2634011543756
## 123  0.529263270 -0.008144690446  0.19215801501  1.2259559 0.2202152756321
## attr(,"call")
## localmoran(x = ssdata$NOR_TXCURR, listw = nb2listw(sscastnb1, 
##     style = "W", zero.policy = TRUE), zero.policy = TRUE, na.action = na.pass)
## attr(,"class")
## [1] "localmoran" "matrix"     "array"     
## attr(,"quadr")
##          mean    median     pysal
## 1   High-High High-High High-High
## 2   High-High High-High High-High
## 3   High-High High-High High-High
## 4   High-High High-High High-High
## 5   High-High High-High High-High
## 6   High-High High-High High-High
## 7    Low-High High-High  Low-High
## 8   High-High High-High High-High
## 9   High-High High-High High-High
## 10  High-High High-High High-High
## 11  High-High High-High High-High
## 12  High-High High-High High-High
## 13   Low-High  Low-High  Low-High
## 14  High-High High-High High-High
## 15   High-Low High-High  High-Low
## 16  High-High High-High High-High
## 17  High-High High-High High-High
## 18   Low-High High-High  Low-High
## 19  High-High High-High High-High
## 20   Low-High  Low-High  Low-High
## 21   Low-High High-High  Low-High
## 22   Low-High High-High  Low-High
## 23  High-High High-High High-High
## 24  High-High High-High High-High
## 25  High-High High-High High-High
## 26  High-High High-High High-High
## 27  High-High High-High High-High
## 28  High-High High-High High-High
## 29  High-High High-High High-High
## 30  High-High High-High High-High
## 31  High-High High-High High-High
## 32    Low-Low   Low-Low   Low-Low
## 33    Low-Low   Low-Low   Low-Low
## 34    Low-Low   Low-Low   Low-Low
## 35    Low-Low   Low-Low   Low-Low
## 36    Low-Low   Low-Low   Low-Low
## 37    Low-Low   Low-Low   Low-Low
## 38    Low-Low   Low-Low   Low-Low
## 39    Low-Low  High-Low   Low-Low
## 40    Low-Low  Low-High   Low-Low
## 41    Low-Low  High-Low   Low-Low
## 42   Low-High  Low-High  Low-High
## 43   High-Low  High-Low  High-Low
## 44    Low-Low  Low-High   Low-Low
## 45    Low-Low   Low-Low   Low-Low
## 46    Low-Low  Low-High   Low-Low
## 47   High-Low  High-Low  High-Low
## 48  High-High High-High High-High
## 49    Low-Low  Low-High   Low-Low
## 50   High-Low  High-Low  High-Low
## 51    Low-Low   Low-Low   Low-Low
## 52    Low-Low  Low-High   Low-Low
## 53    Low-Low   Low-Low   Low-Low
## 54   Low-High  Low-High  Low-High
## 55   High-Low  High-Low  High-Low
## 56   High-Low  High-Low  High-Low
## 57    Low-Low  Low-High   Low-Low
## 58    Low-Low   Low-Low   Low-Low
## 59    Low-Low  Low-High   Low-Low
## 60    Low-Low  High-Low   Low-Low
## 61    Low-Low   Low-Low   Low-Low
## 62    Low-Low  High-Low   Low-Low
## 63    Low-Low   Low-Low   Low-Low
## 64    Low-Low   Low-Low   Low-Low
## 65   High-Low  High-Low  High-Low
## 66    Low-Low   Low-Low   Low-Low
## 67    Low-Low   Low-Low   Low-Low
## 68    Low-Low  High-Low   Low-Low
## 69    Low-Low   Low-Low   Low-Low
## 70    Low-Low   Low-Low   Low-Low
## 71    Low-Low   Low-Low   Low-Low
## 72   High-Low  High-Low  High-Low
## 73    Low-Low   Low-Low   Low-Low
## 74   High-Low  High-Low  High-Low
## 75   High-Low  High-Low  High-Low
## 76    Low-Low   Low-Low   Low-Low
## 77    Low-Low  Low-High   Low-Low
## 78    Low-Low   Low-Low   Low-Low
## 79    Low-Low   Low-Low   Low-Low
## 80    Low-Low  High-Low   Low-Low
## 81    Low-Low  Low-High   Low-Low
## 82   High-Low  High-Low  High-Low
## 83    Low-Low   Low-Low   Low-Low
## 84    Low-Low  High-Low   Low-Low
## 85    Low-Low   Low-Low   Low-Low
## 86    Low-Low  High-Low   Low-Low
## 87    Low-Low   Low-Low   Low-Low
## 88    Low-Low   Low-Low   Low-Low
## 89    Low-Low   Low-Low   Low-Low
## 90    Low-Low   Low-Low   Low-Low
## 91    Low-Low   Low-Low   Low-Low
## 92    Low-Low   Low-Low   Low-Low
## 93    Low-Low   Low-Low   Low-Low
## 94   High-Low  High-Low  High-Low
## 95    Low-Low   Low-Low   Low-Low
## 96    Low-Low   Low-Low   Low-Low
## 97    Low-Low   Low-Low   Low-Low
## 98    Low-Low   Low-Low   Low-Low
## 99    Low-Low   Low-Low   Low-Low
## 100   Low-Low   Low-Low   Low-Low
## 101   Low-Low   Low-Low   Low-Low
## 102   Low-Low  High-Low   Low-Low
## 103   Low-Low  Low-High   Low-Low
## 104   Low-Low   Low-Low   Low-Low
## 105  Low-High  Low-High  Low-High
## 106   Low-Low   Low-Low   Low-Low
## 107  Low-High  Low-High  Low-High
## 108   Low-Low  Low-High   Low-Low
## 109  Low-High  Low-High  Low-High
## 110  Low-High  Low-High  Low-High
## 111  Low-High High-High  Low-High
## 112 High-High High-High High-High
## 113  Low-High High-High  Low-High
## 114 High-High High-High High-High
## 115  High-Low High-High  High-Low
## 116  High-Low  High-Low  High-Low
## 117 High-High High-High High-High
## 118  High-Low High-High  High-Low
## 119   Low-Low  Low-High   Low-Low
## 120  Low-High  Low-High  Low-High
## 121 High-High High-High High-High
## 122 High-High High-High High-High
## 123 High-High High-High High-High
hist(sslmoran)

lmoranmap <- cbind(ssdata, sslmoran)
mdf <- cbind(ssdata, sslmoran, (attributes(sslmoran)$quadr))
tm_shape(lmoranmap) + tm_fill(col = "Ii", style ="quantile",  midpoint = 0, palette = "viridis", title = "Local Moran Observed Statistic") 

tm_shape(lmoranmap) + tm_fill(col = "E.Ii", style ="quantile",  midpoint = 0, palette = "-viridis", title = "Local Moran Expected Statistic") 

tm_shape(lmoranmap) + tm_fill(col = "Z.Ii", style ="quantile",  midpoint = 0, palette = "-viridis", title = "Local Moran Z-Score Statistic")

tm_shape(lmoranmap) + tm_fill(col = "Pr.z....E.Ii..", style ="quantile",  midpoint = 0, palette = "viridis", title = "Local Moran P-Value")

sscastlistw <- nb2listw(sscastnb1, style="W", zero.policy=TRUE)
geary.test(ssdata$NOR_TXCURR, sscastlistw, zero.policy = TRUE)
## 
##  Geary C test under randomisation
## 
## data:  ssdata$NOR_TXCURR 
## weights: sscastlistw 
## 
## Geary C statistic standard deviate = 5.3765, p-value = 0.00000003798
## alternative hypothesis: Expectation greater than statistic
## sample estimates:
## Geary C statistic       Expectation          Variance 
##       0.610042600       1.000000000       0.005260618
geary.mc(ssdata$NOR_TXCURR,sscastlistw, nsim=999, zero.policy = TRUE)
## 
##  Monte-Carlo simulation of Geary C
## 
## data:  ssdata$NOR_TXCURR 
## weights: sscastlistw 
## number of simulations + 1: 1000 
## 
## statistic = 0.61004, observed rank = 1, p-value = 0.001
## alternative hypothesis: greater

Hotspot analysis

sscastlistb <- nb2listw(sscastnb1, style="B", zero.policy=TRUE)
GGis <- globalG.test(ssdata$NOR_TXCURR, sscastlistb, zero.policy = TRUE, alternative = "two.sided")
print(GGis)
## 
##  Getis-Ord global G statistic
## 
## data:  ssdata$NOR_TXCURR 
## weights: sscastlistb 
## 
## standard deviate = 4.9519, p-value = 0.0000007348
## alternative hypothesis: two.sided
## sample estimates:
## Global G statistic        Expectation           Variance 
##      0.06349260630      0.04198320672      0.00001886716
Gis <- localG(ssdata$NOR_TXCURR, sscastlistw, zero.policy = TRUE)
print(Gis)
##   [1]  1.4023330  1.1887789  2.8085119  4.5833048  1.3076790  0.7489160
##   [7]  1.0824950  3.8352450  0.8614194  0.7365262  0.7606685  1.8178513
##  [13]  1.1504692  0.4374518 -0.1489786  1.1988559  4.3644643  2.5609826
##  [19]  2.1973381  0.6673327  2.2698919  1.6221467  2.4917532  2.1773855
##  [25]  3.7127013  1.0360972  4.1859200  1.7907207  1.6201793  5.1650608
##  [31]  1.2121993 -1.6951981 -1.2689635 -0.9406819 -1.1548954 -1.4870953
##  [37] -1.1541951 -1.7103664 -1.3211238 -0.4584639 -1.1656178  2.8064663
##  [43] -0.7105373 -0.2745743 -0.7165964 -0.4339696 -0.7237545  2.4668263
##  [49] -0.2387880 -0.9854741 -0.6654629 -0.1623104 -0.8419595  2.0068068
##  [55] -0.9140990 -1.0794838 -0.3096360 -1.5507657 -0.8581839 -1.8198980
##  [61] -0.6868232 -1.0269352 -1.3299509 -1.5092296 -1.3768758 -0.9386760
##  [67] -1.1023725 -1.0732795 -1.2250267 -1.1561558 -0.8563168 -0.9769950
##  [73] -1.1772107 -1.1002306 -0.9894001 -2.1467160 -0.8309538 -1.0703582
##  [79] -1.0952245 -1.1425662 -0.2426622 -1.0079669 -1.2563033 -1.2243382
##  [85] -1.1988291 -1.6072738 -1.7696402 -1.7693109 -1.2179128 -1.2595835
##  [91] -1.8753429 -1.8053270 -1.1456198 -1.0974466 -1.9708316 -1.3241612
##  [97] -0.8065924 -1.4183492 -1.7017511 -1.9467772 -1.4103553 -0.6997660
## [103] -0.4234433 -1.0082317  0.4714655 -0.6282416  0.7242910 -0.1589480
## [109]  1.9324620  0.3907696  1.2672607  0.5843822  1.0560697  0.2214490
## [115] -0.3263761 -0.8933041  0.6567793 -0.4363135 -0.2224141  0.3733870
## [121]  2.0141640  1.1183884  1.2259559
## attr(,"internals")
##                  Gi       E(Gi)          V(Gi)      Z(Gi)  Pr(z != E(Gi))
##   [1,] 0.0130738227 0.008196721 0.000012095430  1.4023330 0.1608158453402
##   [2,] 0.0131376894 0.008196721 0.000017275154  1.1887789 0.2345266842014
##   [3,] 0.0198717287 0.008196721 0.000017280718  2.8085119 0.0049771040414
##   [4,] 0.0254654129 0.008196721 0.000014195850  4.5833048 0.0000045768422
##   [5,] 0.0131335947 0.008196721 0.000014252849  1.3076790 0.1909822112521
##   [6,] 0.0116851451 0.008196721 0.000021696631  0.7489160 0.4539078643927
##   [7,] 0.0122618151 0.008196721 0.000014102282  1.0824950 0.2790326542547
##   [8,] 0.0224384441 0.008196721 0.000013789193  3.8352450 0.0001254392000
##   [9,] 0.0117686221 0.008196721 0.000017193705  0.8614194 0.3890070965404
##  [10,] 0.0112492970 0.008196721 0.000017177361  0.7365262 0.4614105042275
##  [11,] 0.0132424884 0.008196721 0.000044001104  0.7606685 0.4468551106494
##  [12,] 0.0145077804 0.008196721 0.000012052795  1.8178513 0.0690868588962
##  [13,] 0.0129126950 0.008196721 0.000016803228  1.1504692 0.2499506550092
##  [14,] 0.0105583458 0.008196721 0.000029144813  0.4374518 0.6617837272740
##  [15,] 0.0075028538 0.008196721 0.000021692283 -0.1489786 0.8815705047423
##  [16,] 0.0137832521 0.008196721 0.000021714528  1.1988559 0.2305839685432
##  [17,] 0.0258991127 0.008196721 0.000016451367  4.3644643 0.0000127434661
##  [18,] 0.0164563661 0.008196721 0.000010401823  2.5609826 0.0104376569128
##  [19,] 0.0184020563 0.008196721 0.000021570527  2.1973381 0.0279963073197
##  [20,] 0.0112619816 0.008196721 0.000021098419  0.6673327 0.5045596417553
##  [21,] 0.0155208874 0.008196721 0.000010411325  2.2698919 0.0232141412705
##  [22,] 0.0157138147 0.008196721 0.000021474329  1.6221467 0.1047719322864
##  [23,] 0.0157310831 0.008196721 0.000009142877  2.4917532 0.0127114316274
##  [24,] 0.0164250753 0.008196721 0.000014280890  2.1773855 0.0294518151042
##  [25,] 0.0236013952 0.008196721 0.000017215712  3.7127013 0.0002050588576
##  [26,] 0.0121121752 0.008196721 0.000014281151  1.0360972 0.3001567987920
##  [27,] 0.0277156464 0.008196721 0.000021743525  4.1859200 0.0000284013372
##  [28,] 0.0165185166 0.008196721 0.000021596251  1.7907207 0.0733381308427
##  [29,] 0.0143140912 0.008196721 0.000014256185  1.6201793 0.1051937713384
##  [30,] 0.0296224528 0.008196721 0.000017207605  5.1650608 0.0000002403607
##  [31,] 0.0116833807 0.008196721 0.000008273152  1.2121993 0.2254360931278
##  [32,] 0.0012327248 0.008196721 0.000016876258 -1.6951981 0.0900378429644
##  [33,] 0.0023662738 0.008196721 0.000021110840 -1.2689635 0.2044540765636
##  [34,] 0.0031890616 0.008196721 0.000028338971 -0.9406819 0.3468679126656
##  [35,] 0.0020310182 0.008196721 0.000028502299 -1.1548954 0.2481332545047
##  [36,] 0.0026475162 0.008196721 0.000013924639 -1.4870953 0.1369896382763
##  [37,] 0.0045063055 0.008196721 0.000010223325 -1.1541951 0.2484202099147
##  [38,] 0.0023293255 0.008196721 0.000011768265 -1.7103664 0.0871981361078
##  [39,] 0.0032337174 0.008196721 0.000014112443 -1.3211238 0.1864600830176
##  [40,] 0.0051867937 0.008196721 0.000043102444 -0.4584639 0.6466192282599
##  [41,] 0.0041639600 0.008196721 0.000011969961 -1.1656178 0.2437690159198
##  [42,] 0.0197051854 0.008196721 0.000016815704  2.8064663 0.0050088148951
##  [43,] 0.0018159417 0.008196721 0.000080644318 -0.7105373 0.4773710307771
##  [44,] 0.0067350428 0.008196721 0.000028338971 -0.2745743 0.7836433366700
##  [45,] 0.0048985030 0.008196721 0.000021184080 -0.7165964 0.4736231956126
##  [46,] 0.0064156543 0.008196721 0.000016843871 -0.4339696 0.6643105164214
##  [47,] 0.0043100453 0.008196721 0.000028838598 -0.7237545 0.4692164739028
##  [48,] 0.0184008808 0.008196721 0.000017111076  2.4668263 0.0136316451239
##  [49,] 0.0069255486 0.008196721 0.000028338971 -0.2387880 0.8112699620346
##  [50,] 0.0044812274 0.008196721 0.000014214864 -0.9854741 0.3243912701523
##  [51,] 0.0038226744 0.008196721 0.000043203523 -0.6654629 0.5057544601694
##  [52,] 0.0074482126 0.008196721 0.000021266752 -0.1623104 0.8710614460200
##  [53,] 0.0042997009 0.008196721 0.000021423131 -0.8419595 0.3998106071470
##  [54,] 0.0138956281 0.008196721 0.000008064399  2.0068068 0.0447702424843
##  [55,] 0.0044034363 0.008196721 0.000017220448 -0.9140990 0.3606648056928
##  [56,] 0.0031810998 0.008196721 0.000021588239 -1.0794838 0.2803720900510
##  [57,] 0.0067696966 0.008196721 0.000021240273 -0.3096360 0.7568377846068
##  [58,] 0.0010634376 0.008196721 0.000021158584 -1.5507657 0.1209578337871
##  [59,] 0.0052491254 0.008196721 0.000011797098 -0.8581839 0.3907909232812
##  [60,] 0.0006761049 0.008196721 0.000017077047 -1.8198980 0.0687745340556
##  [61,] 0.0050436403 0.008196721 0.000021075621 -0.6868232 0.4921941593385
##  [62,] 0.0039620445 0.008196721 0.000017004133 -1.0269352 0.3044509877977
##  [63,] 0.0020426300 0.008196721 0.000021411972 -1.3299509 0.1835344562680
##  [64,] 0.0025742454 0.008196721 0.000013878565 -1.5092296 0.1312401097096
##  [65,] 0.0017896203 0.008196721 0.000021653774 -1.3768758 0.1685506406435
##  [66,] 0.0038725467 0.008196721 0.000021221448 -0.9386760 0.3478971116224
##  [67,] 0.0031106763 0.008196721 0.000021286471 -1.1023725 0.2702997645785
##  [68,] 0.0047404194 0.008196721 0.000010370453 -1.0732795 0.2831457464410
##  [69,] 0.0039825165 0.008196721 0.000011834235 -1.2250267 0.2205651680142
##  [70,] 0.0042307753 0.008196721 0.000011766868 -1.1561558 0.2476174385511
##  [71,] 0.0042594861 0.008196721 0.000021140432 -0.8563168 0.3918225963041
##  [72,] 0.0029141293 0.008196721 0.000029235427 -0.9769950 0.3285716145660
##  [73,] 0.0038183277 0.008196721 0.000013833154 -1.1772107 0.2391114269497
##  [74,] 0.0030737402 0.008196721 0.000021680940 -1.1002306 0.2712316804429
##  [75,] 0.0035858303 0.008196721 0.000021718297 -0.9894001 0.3224674024121
##  [76,] 0.0020883682 0.008196721 0.000008096533 -2.1467160 0.0318159016602
##  [77,] 0.0053438992 0.008196721 0.000011786792 -0.8309538 0.4059997323292
##  [78,] 0.0032662067 0.008196721 0.000021219065 -1.0703582 0.2844581187779
##  [79,] 0.0037040044 0.008196721 0.000016827197 -1.0952245 0.2734182764771
##  [80,] 0.0042481448 0.008196721 0.000011943137 -1.1425662 0.2532187625257
##  [81,] 0.0072948707 0.008196721 0.000013812271 -0.2426622 0.8082671268726
##  [82,] 0.0040181164 0.008196721 0.000017185811 -1.0079669 0.3134703310049
##  [83,] 0.0014998196 0.008196721 0.000028415733 -1.2563033 0.2090060386896
##  [84,] 0.0025203742 0.008196721 0.000021494883 -1.2243382 0.2208246627724
##  [85,] 0.0026611627 0.008196721 0.000021321038 -1.1988291 0.2305943995008
##  [86,] 0.0015638975 0.008196721 0.000017030100 -1.6072738 0.1079943173271
##  [87,] 0.0016144767 0.008196721 0.000013834966 -1.7696402 0.0767870976051
##  [88,] 0.0016136113 0.008196721 0.000013843756 -1.7693109 0.0768420109934
##  [89,] 0.0026055040 0.008196721 0.000021075621 -1.2179128 0.2232571261932
##  [90,] 0.0014914094 0.008196721 0.000028338971 -1.2595835 0.2078196617412
##  [91,] 0.0017337789 0.008196721 0.000011876794 -1.8753429 0.0607455715388
##  [92,] 0.0020117888 0.008196721 0.000011737029 -1.8053270 0.0710235207650
##  [93,] 0.0042719028 0.008196721 0.000011737029 -1.1456198 0.2519525106859
##  [94,] 0.0022778257 0.008196721 0.000029088048 -1.0974466 0.2724462122745
##  [95,] 0.0028897615 0.008196721 0.000007250911 -1.9708316 0.0487431412103
##  [96,] 0.0036602316 0.008196721 0.000011737029 -1.3241612 0.1854495224579
##  [97,] 0.0029025136 0.008196721 0.000043081791 -0.8065924 0.4199013806112
##  [98,] 0.0016762825 0.008196721 0.000021134273 -1.4183492 0.1560888325270
##  [99,] 0.0003633390 0.008196721 0.000021188811 -1.7017511 0.0888020455516
## [100,] 0.0015271897 0.008196721 0.000011737029 -1.9467772 0.0515614513102
## [101,] 0.0033555751 0.008196721 0.000011782551 -1.4103553 0.1584348038729
## [102,] 0.0049479007 0.008196721 0.000021554890 -0.6997660 0.4840734551855
## [103,] 0.0064601012 0.008196721 0.000016819733 -0.4234433 0.6719718542450
## [104,] 0.0028282658 0.008196721 0.000028351628 -1.0082317 0.3133432364167
## [105,] 0.0103611343 0.008196721 0.000021075621  0.4714655 0.6373083277561
## [106,] 0.0040799227 0.008196721 0.000042940382 -0.6282416 0.5298457064391
## [107,] 0.0109014931 0.008196721 0.000013945530  0.7242910 0.4688870768594
## [108,] 0.0076859990 0.008196721 0.000010324279 -0.1589480 0.8737098788773
## [109,] 0.0154332467 0.008196721 0.000014022914  1.9324620 0.0533025136999
## [110,] 0.0095433779 0.008196721 0.000011876051  0.3907696 0.6959675118241
## [111,] 0.0140673980 0.008196721 0.000021460724  1.2672607 0.2050620754513
## [112,] 0.0106258057 0.008196721 0.000017277914  0.5843822 0.5589632305307
## [113,] 0.0138807140 0.008196721 0.000028968221  1.0560697 0.2909363503425
## [114,] 0.0088708079 0.008196721 0.000009265819  0.2214490 0.8247428182119
## [115,] 0.0070690985 0.008196721 0.000011936886 -0.3263761 0.7441398363457
## [116,] 0.0048261719 0.008196721 0.000014236485 -0.8933041 0.3716943550869
## [117,] 0.0109175072 0.008196721 0.000017161286  0.6567793 0.5113228454889
## [118,] 0.0063865137 0.008196721 0.000017213113 -0.4363135 0.6626092724731
## [119,] 0.0061135242 0.008196721 0.000087727547 -0.2224141 0.8239915098185
## [120,] 0.0099108733 0.008196721 0.000021075621  0.3733870 0.7088604661011
## [121,] 0.0165514569 0.008196721 0.000017205836  2.0141640 0.0439923206612
## [122,] 0.0133945617 0.008196721 0.000021600340  1.1183884 0.2634011543756
## [123,] 0.0132932269 0.008196721 0.000017282052  1.2259559 0.2202152756321
## attr(,"cluster")
##   [1] High High High High High High Low  High High High High High Low  High High
##  [16] High High Low  High Low  Low  Low  High High High High High High High High
##  [31] High Low  Low  Low  Low  Low  Low  Low  Low  Low  Low  Low  High Low  Low 
##  [46] Low  High High Low  High Low  Low  Low  Low  High High Low  Low  Low  Low 
##  [61] Low  Low  Low  Low  High Low  Low  Low  Low  Low  Low  High Low  High High
##  [76] Low  Low  Low  Low  Low  Low  High Low  Low  Low  Low  Low  Low  Low  Low 
##  [91] Low  Low  Low  High Low  Low  Low  Low  Low  Low  Low  Low  Low  Low  Low 
## [106] Low  Low  Low  Low  Low  Low  High Low  High High High High High Low  Low 
## [121] High High High
## Levels: Low High
## attr(,"gstari")
## [1] FALSE
## attr(,"call")
## localG(x = ssdata$NOR_TXCURR, listw = sscastlistw, zero.policy = TRUE)
## attr(,"class")
## [1] "localG"
Gimap <- cbind(ssdata, as.matrix(Gis), as.matrix(attributes(Gis)$internals), as.matrix(attributes(Gis)$cluster))
Gimap <- Gimap %>% rename('prz'= 'Pr.z....E.Gi..')
Gimap <- Gimap %>% rename_at('as.matrix.attributes.Gis..cluster.', ~'cluster')
tm_shape(Gimap) + tm_fill(col = "as.matrix.Gis.", palette = "-viridis", midpoint = NA, title = "G*") + tm_legend(position = c("centre", "top"))

tm_shape(Gimap) + tm_fill(col = "Gi", palette = "-viridis", title = "Local G Observed Statistic")

tm_shape(Gimap) + tm_fill(col = "prz", style ="quantile", palette = "-viridis", title = "Local G P-Value") + tm_layout(legend.position = c(0.45, 0.65))

Aspatial multivariate analysis using Poisson distribution

Poisson regression

poisson_output <-glm(formula = NOR_TXCURR ~ POP_DEN + PHF_COV + AHF_COV + POP_POV + HEALTH_DEP + EDUC_DEP + WORK_DEP + SEC_SHOCK,
data = rschdata, family = poisson)
## Warning in dpois(y, mu, log = TRUE): non-integer x = 8481.401340
## Warning in dpois(y, mu, log = TRUE): non-integer x = 5871.310184
## Warning in dpois(y, mu, log = TRUE): non-integer x = 4959.670205
## Warning in dpois(y, mu, log = TRUE): non-integer x = 5471.895840
## Warning in dpois(y, mu, log = TRUE): non-integer x = 8420.317920
## Warning in dpois(y, mu, log = TRUE): non-integer x = 2038.665138
## Warning in dpois(y, mu, log = TRUE): non-integer x = 7014.040846
## Warning in dpois(y, mu, log = TRUE): non-integer x = 10895.600778
## Warning in dpois(y, mu, log = TRUE): non-integer x = 2953.330215
## Warning in dpois(y, mu, log = TRUE): non-integer x = 2460.296866
## Warning in dpois(y, mu, log = TRUE): non-integer x = 4479.285188
## Warning in dpois(y, mu, log = TRUE): non-integer x = 3502.882481
## Warning in dpois(y, mu, log = TRUE): non-integer x = 2691.169494
## Warning in dpois(y, mu, log = TRUE): non-integer x = 4673.598586
## Warning in dpois(y, mu, log = TRUE): non-integer x = 10005.176030
## Warning in dpois(y, mu, log = TRUE): non-integer x = 8248.138284
## Warning in dpois(y, mu, log = TRUE): non-integer x = 3740.698556
## Warning in dpois(y, mu, log = TRUE): non-integer x = 5035.542691
## Warning in dpois(y, mu, log = TRUE): non-integer x = 3817.352694
## Warning in dpois(y, mu, log = TRUE): non-integer x = 3264.867055
## Warning in dpois(y, mu, log = TRUE): non-integer x = 6173.078593
## Warning in dpois(y, mu, log = TRUE): non-integer x = 6754.753217
## Warning in dpois(y, mu, log = TRUE): non-integer x = 3358.161900
## Warning in dpois(y, mu, log = TRUE): non-integer x = 2865.751974
## Warning in dpois(y, mu, log = TRUE): non-integer x = 3451.472782
## Warning in dpois(y, mu, log = TRUE): non-integer x = 1818.172777
## Warning in dpois(y, mu, log = TRUE): non-integer x = 2166.468698
## Warning in dpois(y, mu, log = TRUE): non-integer x = 1204.255817
## Warning in dpois(y, mu, log = TRUE): non-integer x = 7128.989005
## Warning in dpois(y, mu, log = TRUE): non-integer x = 3038.513644
## Warning in dpois(y, mu, log = TRUE): non-integer x = 3410.668567
## Warning in dpois(y, mu, log = TRUE): non-integer x = 13298.755187
## Warning in dpois(y, mu, log = TRUE): non-integer x = 2854.128317
## Warning in dpois(y, mu, log = TRUE): non-integer x = 3165.654120
## Warning in dpois(y, mu, log = TRUE): non-integer x = 2378.993253
## Warning in dpois(y, mu, log = TRUE): non-integer x = 4920.580700
## Warning in dpois(y, mu, log = TRUE): non-integer x = 3083.709932
## Warning in dpois(y, mu, log = TRUE): non-integer x = 4380.800637
## Warning in dpois(y, mu, log = TRUE): non-integer x = 1948.638019
## Warning in dpois(y, mu, log = TRUE): non-integer x = 3461.228178
## Warning in dpois(y, mu, log = TRUE): non-integer x = 2038.418543
## Warning in dpois(y, mu, log = TRUE): non-integer x = 989.605252
## Warning in dpois(y, mu, log = TRUE): non-integer x = 2097.851365
## Warning in dpois(y, mu, log = TRUE): non-integer x = 1317.327670
## Warning in dpois(y, mu, log = TRUE): non-integer x = 1635.274435
## Warning in dpois(y, mu, log = TRUE): non-integer x = 1512.785015
## Warning in dpois(y, mu, log = TRUE): non-integer x = 3176.829154
## Warning in dpois(y, mu, log = TRUE): non-integer x = 6305.475941
## Warning in dpois(y, mu, log = TRUE): non-integer x = 3197.854925
## Warning in dpois(y, mu, log = TRUE): non-integer x = 2766.486275
## Warning in dpois(y, mu, log = TRUE): non-integer x = 3229.232904
## Warning in dpois(y, mu, log = TRUE): non-integer x = 1389.514188
## Warning in dpois(y, mu, log = TRUE): non-integer x = 788.023557
## Warning in dpois(y, mu, log = TRUE): non-integer x = 2271.334457
## Warning in dpois(y, mu, log = TRUE): non-integer x = 2156.529404
## Warning in dpois(y, mu, log = TRUE): non-integer x = 993.581464
## Warning in dpois(y, mu, log = TRUE): non-integer x = 1706.962657
## Warning in dpois(y, mu, log = TRUE): non-integer x = 2432.547536
## Warning in dpois(y, mu, log = TRUE): non-integer x = 1444.338618
## Warning in dpois(y, mu, log = TRUE): non-integer x = 3187.785293
## Warning in dpois(y, mu, log = TRUE): non-integer x = 2264.150943
## Warning in dpois(y, mu, log = TRUE): non-integer x = 1812.692485
## Warning in dpois(y, mu, log = TRUE): non-integer x = 463.902066
## Warning in dpois(y, mu, log = TRUE): non-integer x = 1557.858493
## Warning in dpois(y, mu, log = TRUE): non-integer x = 2708.026111
## Warning in dpois(y, mu, log = TRUE): non-integer x = 819.818382
## Warning in dpois(y, mu, log = TRUE): non-integer x = 2699.916909
## Warning in dpois(y, mu, log = TRUE): non-integer x = 1681.717813
## Warning in dpois(y, mu, log = TRUE): non-integer x = 483.908382
## Warning in dpois(y, mu, log = TRUE): non-integer x = 987.692308
## Warning in dpois(y, mu, log = TRUE): non-integer x = 2010.017072
## Warning in dpois(y, mu, log = TRUE): non-integer x = 1676.891283
## Warning in dpois(y, mu, log = TRUE): non-integer x = 630.986256
## Warning in dpois(y, mu, log = TRUE): non-integer x = 1274.897249
## Warning in dpois(y, mu, log = TRUE): non-integer x = 627.822492
## Warning in dpois(y, mu, log = TRUE): non-integer x = 301.822941
## Warning in dpois(y, mu, log = TRUE): non-integer x = 521.197867
## Warning in dpois(y, mu, log = TRUE): non-integer x = 736.868232
## Warning in dpois(y, mu, log = TRUE): non-integer x = 969.112947
## Warning in dpois(y, mu, log = TRUE): non-integer x = 774.004910
## Warning in dpois(y, mu, log = TRUE): non-integer x = 633.764206
## Warning in dpois(y, mu, log = TRUE): non-integer x = 668.274128
## Warning in dpois(y, mu, log = TRUE): non-integer x = 605.577233
## Warning in dpois(y, mu, log = TRUE): non-integer x = 655.602657
## Warning in dpois(y, mu, log = TRUE): non-integer x = 402.460159
## Warning in dpois(y, mu, log = TRUE): non-integer x = 435.482057
## Warning in dpois(y, mu, log = TRUE): non-integer x = 208.910468
## Warning in dpois(y, mu, log = TRUE): non-integer x = 404.490875
## Warning in dpois(y, mu, log = TRUE): non-integer x = 330.872183
## Warning in dpois(y, mu, log = TRUE): non-integer x = 456.819556
## Warning in dpois(y, mu, log = TRUE): non-integer x = 551.828945
## Warning in dpois(y, mu, log = TRUE): non-integer x = 376.207687
## Warning in dpois(y, mu, log = TRUE): non-integer x = 327.422414
## Warning in dpois(y, mu, log = TRUE): non-integer x = 542.296508
## Warning in dpois(y, mu, log = TRUE): non-integer x = 395.957829
## Warning in dpois(y, mu, log = TRUE): non-integer x = 205.185600
## Warning in dpois(y, mu, log = TRUE): non-integer x = 402.178212
## Warning in dpois(y, mu, log = TRUE): non-integer x = 127.888291
## Warning in dpois(y, mu, log = TRUE): non-integer x = 455.271689
## Warning in dpois(y, mu, log = TRUE): non-integer x = 237.828891
## Warning in dpois(y, mu, log = TRUE): non-integer x = 422.899080
## Warning in dpois(y, mu, log = TRUE): non-integer x = 214.751508
## Warning in dpois(y, mu, log = TRUE): non-integer x = 306.464857
## Warning in dpois(y, mu, log = TRUE): non-integer x = 195.831658
## Warning in dpois(y, mu, log = TRUE): non-integer x = 125.718017
## Warning in dpois(y, mu, log = TRUE): non-integer x = 133.459532
## Warning in dpois(y, mu, log = TRUE): non-integer x = 175.237327
## Warning in dpois(y, mu, log = TRUE): non-integer x = 82.428390
## Warning in dpois(y, mu, log = TRUE): non-integer x = 115.570751
## Warning in dpois(y, mu, log = TRUE): non-integer x = 33.879776
print(summary(poisson_output))
## 
## Call:
## glm(formula = NOR_TXCURR ~ POP_DEN + PHF_COV + AHF_COV + POP_POV + 
##     HEALTH_DEP + EDUC_DEP + WORK_DEP + SEC_SHOCK, family = poisson, 
##     data = rschdata)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -97.790  -17.152   -0.886   17.732   83.428  
## 
## Coefficients:
##                   Estimate     Std. Error z value            Pr(>|z|)    
## (Intercept)  8.22781639879  0.00628298986 1309.54 <0.0000000000000002 ***
## POP_DEN      0.00006018465  0.00000219290   27.45 <0.0000000000000002 ***
## PHF_COV      0.01568730882  0.00016649608   94.22 <0.0000000000000002 ***
## AHF_COV     -0.00618443118  0.00003826263 -161.63 <0.0000000000000002 ***
## POP_POV     -0.00000317098  0.00000005951  -53.28 <0.0000000000000002 ***
## HEALTH_DEP   0.00000701215  0.00000018534   37.83 <0.0000000000000002 ***
## EDUC_DEP    -0.00003776223  0.00000031499 -119.89 <0.0000000000000002 ***
## WORK_DEP     0.00002833237  0.00000026063  108.71 <0.0000000000000002 ***
## SEC_SHOCK   -0.00006747833  0.00000042840 -157.51 <0.0000000000000002 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 304081  on 122  degrees of freedom
## Residual deviance: 136424  on 114  degrees of freedom
## AIC: Inf
## 
## Number of Fisher Scoring iterations: 5
vif_values <- vif(poisson_output)
vif_values
##    POP_DEN    PHF_COV    AHF_COV    POP_POV HEALTH_DEP   EDUC_DEP   WORK_DEP 
##   1.614147   1.360960   1.072502   8.350692  22.275535   3.684154  14.762289 
##  SEC_SHOCK 
##   3.907127
barplot(vif_values, main = "VIF Values", horiz = TRUE, col = "black", cex.names = 0.5)+ abline(v = 5, lwd = 3, lty = 2)

## numeric(0)
plot(poisson_output)

poisson_output1 <-glm(formula = NOR_TXCURR ~ POP_DEN + PHF_COV + AHF_COV + POP_POV + EDUC_DEP + WORK_DEP + SEC_SHOCK,
data = rschdata, family = poisson)
## Warning in dpois(y, mu, log = TRUE): non-integer x = 8481.401340
## Warning in dpois(y, mu, log = TRUE): non-integer x = 5871.310184
## Warning in dpois(y, mu, log = TRUE): non-integer x = 4959.670205
## Warning in dpois(y, mu, log = TRUE): non-integer x = 5471.895840
## Warning in dpois(y, mu, log = TRUE): non-integer x = 8420.317920
## Warning in dpois(y, mu, log = TRUE): non-integer x = 2038.665138
## Warning in dpois(y, mu, log = TRUE): non-integer x = 7014.040846
## Warning in dpois(y, mu, log = TRUE): non-integer x = 10895.600778
## Warning in dpois(y, mu, log = TRUE): non-integer x = 2953.330215
## Warning in dpois(y, mu, log = TRUE): non-integer x = 2460.296866
## Warning in dpois(y, mu, log = TRUE): non-integer x = 4479.285188
## Warning in dpois(y, mu, log = TRUE): non-integer x = 3502.882481
## Warning in dpois(y, mu, log = TRUE): non-integer x = 2691.169494
## Warning in dpois(y, mu, log = TRUE): non-integer x = 4673.598586
## Warning in dpois(y, mu, log = TRUE): non-integer x = 10005.176030
## Warning in dpois(y, mu, log = TRUE): non-integer x = 8248.138284
## Warning in dpois(y, mu, log = TRUE): non-integer x = 3740.698556
## Warning in dpois(y, mu, log = TRUE): non-integer x = 5035.542691
## Warning in dpois(y, mu, log = TRUE): non-integer x = 3817.352694
## Warning in dpois(y, mu, log = TRUE): non-integer x = 3264.867055
## Warning in dpois(y, mu, log = TRUE): non-integer x = 6173.078593
## Warning in dpois(y, mu, log = TRUE): non-integer x = 6754.753217
## Warning in dpois(y, mu, log = TRUE): non-integer x = 3358.161900
## Warning in dpois(y, mu, log = TRUE): non-integer x = 2865.751974
## Warning in dpois(y, mu, log = TRUE): non-integer x = 3451.472782
## Warning in dpois(y, mu, log = TRUE): non-integer x = 1818.172777
## Warning in dpois(y, mu, log = TRUE): non-integer x = 2166.468698
## Warning in dpois(y, mu, log = TRUE): non-integer x = 1204.255817
## Warning in dpois(y, mu, log = TRUE): non-integer x = 7128.989005
## Warning in dpois(y, mu, log = TRUE): non-integer x = 3038.513644
## Warning in dpois(y, mu, log = TRUE): non-integer x = 3410.668567
## Warning in dpois(y, mu, log = TRUE): non-integer x = 13298.755187
## Warning in dpois(y, mu, log = TRUE): non-integer x = 2854.128317
## Warning in dpois(y, mu, log = TRUE): non-integer x = 3165.654120
## Warning in dpois(y, mu, log = TRUE): non-integer x = 2378.993253
## Warning in dpois(y, mu, log = TRUE): non-integer x = 4920.580700
## Warning in dpois(y, mu, log = TRUE): non-integer x = 3083.709932
## Warning in dpois(y, mu, log = TRUE): non-integer x = 4380.800637
## Warning in dpois(y, mu, log = TRUE): non-integer x = 1948.638019
## Warning in dpois(y, mu, log = TRUE): non-integer x = 3461.228178
## Warning in dpois(y, mu, log = TRUE): non-integer x = 2038.418543
## Warning in dpois(y, mu, log = TRUE): non-integer x = 989.605252
## Warning in dpois(y, mu, log = TRUE): non-integer x = 2097.851365
## Warning in dpois(y, mu, log = TRUE): non-integer x = 1317.327670
## Warning in dpois(y, mu, log = TRUE): non-integer x = 1635.274435
## Warning in dpois(y, mu, log = TRUE): non-integer x = 1512.785015
## Warning in dpois(y, mu, log = TRUE): non-integer x = 3176.829154
## Warning in dpois(y, mu, log = TRUE): non-integer x = 6305.475941
## Warning in dpois(y, mu, log = TRUE): non-integer x = 3197.854925
## Warning in dpois(y, mu, log = TRUE): non-integer x = 2766.486275
## Warning in dpois(y, mu, log = TRUE): non-integer x = 3229.232904
## Warning in dpois(y, mu, log = TRUE): non-integer x = 1389.514188
## Warning in dpois(y, mu, log = TRUE): non-integer x = 788.023557
## Warning in dpois(y, mu, log = TRUE): non-integer x = 2271.334457
## Warning in dpois(y, mu, log = TRUE): non-integer x = 2156.529404
## Warning in dpois(y, mu, log = TRUE): non-integer x = 993.581464
## Warning in dpois(y, mu, log = TRUE): non-integer x = 1706.962657
## Warning in dpois(y, mu, log = TRUE): non-integer x = 2432.547536
## Warning in dpois(y, mu, log = TRUE): non-integer x = 1444.338618
## Warning in dpois(y, mu, log = TRUE): non-integer x = 3187.785293
## Warning in dpois(y, mu, log = TRUE): non-integer x = 2264.150943
## Warning in dpois(y, mu, log = TRUE): non-integer x = 1812.692485
## Warning in dpois(y, mu, log = TRUE): non-integer x = 463.902066
## Warning in dpois(y, mu, log = TRUE): non-integer x = 1557.858493
## Warning in dpois(y, mu, log = TRUE): non-integer x = 2708.026111
## Warning in dpois(y, mu, log = TRUE): non-integer x = 819.818382
## Warning in dpois(y, mu, log = TRUE): non-integer x = 2699.916909
## Warning in dpois(y, mu, log = TRUE): non-integer x = 1681.717813
## Warning in dpois(y, mu, log = TRUE): non-integer x = 483.908382
## Warning in dpois(y, mu, log = TRUE): non-integer x = 987.692308
## Warning in dpois(y, mu, log = TRUE): non-integer x = 2010.017072
## Warning in dpois(y, mu, log = TRUE): non-integer x = 1676.891283
## Warning in dpois(y, mu, log = TRUE): non-integer x = 630.986256
## Warning in dpois(y, mu, log = TRUE): non-integer x = 1274.897249
## Warning in dpois(y, mu, log = TRUE): non-integer x = 627.822492
## Warning in dpois(y, mu, log = TRUE): non-integer x = 301.822941
## Warning in dpois(y, mu, log = TRUE): non-integer x = 521.197867
## Warning in dpois(y, mu, log = TRUE): non-integer x = 736.868232
## Warning in dpois(y, mu, log = TRUE): non-integer x = 969.112947
## Warning in dpois(y, mu, log = TRUE): non-integer x = 774.004910
## Warning in dpois(y, mu, log = TRUE): non-integer x = 633.764206
## Warning in dpois(y, mu, log = TRUE): non-integer x = 668.274128
## Warning in dpois(y, mu, log = TRUE): non-integer x = 605.577233
## Warning in dpois(y, mu, log = TRUE): non-integer x = 655.602657
## Warning in dpois(y, mu, log = TRUE): non-integer x = 402.460159
## Warning in dpois(y, mu, log = TRUE): non-integer x = 435.482057
## Warning in dpois(y, mu, log = TRUE): non-integer x = 208.910468
## Warning in dpois(y, mu, log = TRUE): non-integer x = 404.490875
## Warning in dpois(y, mu, log = TRUE): non-integer x = 330.872183
## Warning in dpois(y, mu, log = TRUE): non-integer x = 456.819556
## Warning in dpois(y, mu, log = TRUE): non-integer x = 551.828945
## Warning in dpois(y, mu, log = TRUE): non-integer x = 376.207687
## Warning in dpois(y, mu, log = TRUE): non-integer x = 327.422414
## Warning in dpois(y, mu, log = TRUE): non-integer x = 542.296508
## Warning in dpois(y, mu, log = TRUE): non-integer x = 395.957829
## Warning in dpois(y, mu, log = TRUE): non-integer x = 205.185600
## Warning in dpois(y, mu, log = TRUE): non-integer x = 402.178212
## Warning in dpois(y, mu, log = TRUE): non-integer x = 127.888291
## Warning in dpois(y, mu, log = TRUE): non-integer x = 455.271689
## Warning in dpois(y, mu, log = TRUE): non-integer x = 237.828891
## Warning in dpois(y, mu, log = TRUE): non-integer x = 422.899080
## Warning in dpois(y, mu, log = TRUE): non-integer x = 214.751508
## Warning in dpois(y, mu, log = TRUE): non-integer x = 306.464857
## Warning in dpois(y, mu, log = TRUE): non-integer x = 195.831658
## Warning in dpois(y, mu, log = TRUE): non-integer x = 125.718017
## Warning in dpois(y, mu, log = TRUE): non-integer x = 133.459532
## Warning in dpois(y, mu, log = TRUE): non-integer x = 175.237327
## Warning in dpois(y, mu, log = TRUE): non-integer x = 82.428390
## Warning in dpois(y, mu, log = TRUE): non-integer x = 115.570751
## Warning in dpois(y, mu, log = TRUE): non-integer x = 33.879776
print(summary(poisson_output1))
## 
## Call:
## glm(formula = NOR_TXCURR ~ POP_DEN + PHF_COV + AHF_COV + POP_POV + 
##     EDUC_DEP + WORK_DEP + SEC_SHOCK, family = poisson, data = rschdata)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -98.621  -19.121    0.351   18.421   81.594  
## 
## Coefficients:
##                   Estimate     Std. Error z value            Pr(>|z|)    
## (Intercept)  8.31858168385  0.00583431141 1425.80 <0.0000000000000002 ***
## POP_DEN      0.00008990058  0.00000203674   44.14 <0.0000000000000002 ***
## PHF_COV      0.01400305796  0.00016128267   86.82 <0.0000000000000002 ***
## AHF_COV     -0.00626078615  0.00003784892 -165.41 <0.0000000000000002 ***
## POP_POV     -0.00000216883  0.00000005202  -41.69 <0.0000000000000002 ***
## EDUC_DEP    -0.00003102203  0.00000025882 -119.86 <0.0000000000000002 ***
## WORK_DEP     0.00003336394  0.00000023380  142.70 <0.0000000000000002 ***
## SEC_SHOCK   -0.00006514925  0.00000042239 -154.24 <0.0000000000000002 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 304081  on 122  degrees of freedom
## Residual deviance: 137907  on 115  degrees of freedom
## AIC: Inf
## 
## Number of Fisher Scoring iterations: 5
vif_values1 <- vif(poisson_output1)
vif_values1
##   POP_DEN   PHF_COV   AHF_COV   POP_POV  EDUC_DEP  WORK_DEP SEC_SHOCK 
##  1.376051  1.254362  1.062385  6.972575  2.511460 12.106627  3.837199
library(MASS)
## Warning: package 'MASS' was built under R version 4.2.3
## 
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
## 
##     select
gnb <- glm(NOR_TXCURR ~ POP_DEN + PHF_COV + AHF_COV + POP_POV + HEALTH_DEP + EDUC_DEP + WORK_DEP + SEC_SHOCK, family = negative.binomial(2), data = rschdata)
print(summary(gnb))
## 
## Call:
## glm(formula = NOR_TXCURR ~ POP_DEN + PHF_COV + AHF_COV + POP_POV + 
##     HEALTH_DEP + EDUC_DEP + WORK_DEP + SEC_SHOCK, family = negative.binomial(2), 
##     data = rschdata)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -5.5655  -1.0226  -0.1408   0.5074   2.2560  
## 
## Coefficients:
##                 Estimate   Std. Error t value             Pr(>|t|)    
## (Intercept)  8.055704222  0.248862370  32.370 < 0.0000000000000002 ***
## POP_DEN      0.000213907  0.000099261   2.155              0.03327 *  
## PHF_COV      0.006443845  0.007806149   0.825              0.41082    
## AHF_COV     -0.003067233  0.000383892  -7.990      0.0000000000012 ***
## POP_POV     -0.000002114  0.000001927  -1.097              0.27508    
## HEALTH_DEP   0.000002879  0.000006242   0.461              0.64544    
## EDUC_DEP    -0.000027998  0.000010417  -2.688              0.00827 ** 
## WORK_DEP     0.000018453  0.000006549   2.818              0.00570 ** 
## SEC_SHOCK   -0.000040491  0.000012146  -3.334              0.00116 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(2) family taken to be 1.286706)
## 
##     Null deviance: 628.15  on 122  degrees of freedom
## Residual deviance: 466.13  on 114  degrees of freedom
## AIC: 2204.7
## 
## Number of Fisher Scoring iterations: 24
pchisq(466.25, 114, lower.tail = FALSE)
## [1] 0.00000000000000000000000000000000000000000004043444
lm.morantest(gnb, sscastlistw, alternative = "two.sided")
## 
##  Global Moran I for regression residuals
## 
## data:  
## model: glm(formula = NOR_TXCURR ~ POP_DEN + PHF_COV + AHF_COV + POP_POV
## + HEALTH_DEP + EDUC_DEP + WORK_DEP + SEC_SHOCK, family =
## negative.binomial(2), data = rschdata)
## weights: sscastlistw
## 
## Moran I statistic standard deviate = 9.6301, p-value <
## 0.00000000000000022
## alternative hypothesis: two.sided
## sample estimates:
## Observed Moran I      Expectation         Variance 
##      0.559071077     -0.008799756      0.003477290
lm.LMtests(gnb, sscastlistw, test = c("LMerr", "LMlag", "RLMerr", "RLMlag", "SARMA"))
## 
##  Lagrange multiplier diagnostics for spatial dependence
## 
## data:  
## model: glm(formula = NOR_TXCURR ~ POP_DEN + PHF_COV + AHF_COV + POP_POV
## + HEALTH_DEP + EDUC_DEP + WORK_DEP + SEC_SHOCK, family =
## negative.binomial(2), data = rschdata)
## weights: sscastlistw
## 
## LMerr = 90.124, df = 1, p-value < 0.00000000000000022
## 
## 
##  Lagrange multiplier diagnostics for spatial dependence
## 
## data:  
## model: glm(formula = NOR_TXCURR ~ POP_DEN + PHF_COV + AHF_COV + POP_POV
## + HEALTH_DEP + EDUC_DEP + WORK_DEP + SEC_SHOCK, family =
## negative.binomial(2), data = rschdata)
## weights: sscastlistw
## 
## LMlag = 2.9942, df = 1, p-value = 0.08356
## 
## 
##  Lagrange multiplier diagnostics for spatial dependence
## 
## data:  
## model: glm(formula = NOR_TXCURR ~ POP_DEN + PHF_COV + AHF_COV + POP_POV
## + HEALTH_DEP + EDUC_DEP + WORK_DEP + SEC_SHOCK, family =
## negative.binomial(2), data = rschdata)
## weights: sscastlistw
## 
## RLMerr = 90.099, df = 1, p-value < 0.00000000000000022
## 
## 
##  Lagrange multiplier diagnostics for spatial dependence
## 
## data:  
## model: glm(formula = NOR_TXCURR ~ POP_DEN + PHF_COV + AHF_COV + POP_POV
## + HEALTH_DEP + EDUC_DEP + WORK_DEP + SEC_SHOCK, family =
## negative.binomial(2), data = rschdata)
## weights: sscastlistw
## 
## RLMlag = 2.9692, df = 1, p-value = 0.08487
## 
## 
##  Lagrange multiplier diagnostics for spatial dependence
## 
## data:  
## model: glm(formula = NOR_TXCURR ~ POP_DEN + PHF_COV + AHF_COV + POP_POV
## + HEALTH_DEP + EDUC_DEP + WORK_DEP + SEC_SHOCK, family =
## negative.binomial(2), data = rschdata)
## weights: sscastlistw
## 
## SARMA = 93.093, df = 2, p-value < 0.00000000000000022

Geographical weighted poisson regression

ssdataspdf <- sf:::as_Spatial(ssdata)
DM<-gw.dist(dp.locat=coordinates(ssdataspdf))
ssdataspdf$NOR_TXCURR <- round(ssdataspdf$NOR_TXCURR, digits = 0)
bw <- bw.ggwr(NOR_TXCURR ~ POP_DEN + PHF_COV + AHF_COV + POP_POV + HEALTH_DEP + EDUC_DEP + WORK_DEP + SEC_SHOCK, data = ssdataspdf, family = "poisson", approach = "AICc", kernel = "exponential", adaptive = TRUE, dMat = DM )
##  Iteration    Log-Likelihood(With bandwidth:  83 )
## =========================
##        0      -2.325e+05 
##        1      -2.779e+09 
##        2      -1.023e+09 
##        3      -3.762e+08 
##        4      -1.385e+08 
##        5      -5.1e+07 
##        6      -1.882e+07 
##        7      -6.974e+06 
##        8      -2.603e+06 
##        9      -9.892e+05 
##       10      -3.912e+05 
##       11      -1.701e+05 
##       12      -9.073e+04 
##       13      -6.528e+04 
##       14      -5.9e+04 
##       15      -5.761e+04 
##       16      -5.728e+04 
##       17      -5.727e+04 
##       18      -5.727e+04 
##       19      -5.727e+04 
## Adaptive bandwidth (number of nearest neighbours): 83 AICc value: 114566.1 
##  Iteration    Log-Likelihood(With bandwidth:  59 )
## =========================
##        0      -1.966e+05 
##        1      -2.149e+07 
##        2      -7.924e+06 
##        3      -2.915e+06 
##        4      -1.082e+06 
##        5      -4.158e+05 
##        6      -1.759e+05 
##        7      -9.096e+04 
##        8      -6.338e+04 
##        9      -5.613e+04 
##       10      -5.437e+04 
##       11      -5.397e+04 
##       12      -5.394e+04 
##       13      -5.394e+04 
##       14      -5.394e+04 
## Adaptive bandwidth (number of nearest neighbours): 59 AICc value: 107922.9 
##  Iteration    Log-Likelihood(With bandwidth:  43 )
## =========================
##        0      -1.612e+05 
##        1      -1.572e+06 
##        2      -5.604e+05 
##        3      -2.083e+05 
##        4      -9.491e+04 
##        5      -6.116e+04 
##        6      -5.173e+04 
##        7      -4.983e+04 
##        8      -4.953e+04 
##        9      -4.949e+04 
##       10      -4.949e+04 
## Adaptive bandwidth (number of nearest neighbours): 43 AICc value: 99024.94 
##  Iteration    Log-Likelihood(With bandwidth:  34 )
## =========================
##        0      -1.412e+05 
##        1      -1.484e+06 
##        2      -5.596e+05 
##        3      -2.169e+05 
##        4      -1.004e+05 
##        5      -6.217e+04 
##        6      -5.001e+04 
##        7      -4.689e+04 
##        8      -4.631e+04 
##        9      -4.619e+04 
##       10      -4.618e+04 
##       11      -4.618e+04 
## Adaptive bandwidth (number of nearest neighbours): 34 AICc value: 92424.72 
##  Iteration    Log-Likelihood(With bandwidth:  27 )
## =========================
##        0      -1.277e+05 
##        1      -1.953e+06 
##        2      -7.668e+05 
##        3      -2.979e+05 
##        4      -1.295e+05 
##        5      -7.1e+04 
##        6      -5.099e+04 
##        7      -4.497e+04 
##        8      -4.367e+04 
##        9      -4.34e+04 
##       10      -4.336e+04 
##       11      -4.336e+04 
## Adaptive bandwidth (number of nearest neighbours): 27 AICc value: 86784.45 
##  Iteration    Log-Likelihood(With bandwidth:  24 )
## =========================
##        0      -1.217e+05 
##        1      -2.497e+06 
##        2      -1.003e+06 
##        3      -3.909e+05 
##        4      -1.653e+05 
##        5      -8.386e+04 
##        6      -5.435e+04 
##        7      -4.475e+04 
##        8      -4.239e+04 
##        9      -4.191e+04 
##       10      -4.181e+04 
##       11      -4.181e+04 
##       12      -4.181e+04 
##       13      -4.181e+04 
## Adaptive bandwidth (number of nearest neighbours): 24 AICc value: 83690.03 
##  Iteration    Log-Likelihood(With bandwidth:  21 )
## =========================
##        0      -1.15e+05 
##        1      -2.547e+06 
##        2      -1.016e+06 
##        3      -3.939e+05 
##        4      -1.658e+05 
##        5      -8.316e+04 
##        6      -5.345e+04 
##        7      -4.363e+04 
##        8      -4.103e+04 
##        9      -4.047e+04 
##       10      -4.035e+04 
##       11      -4.035e+04 
##       12      -4.035e+04 
##       13      -4.035e+04 
## Adaptive bandwidth (number of nearest neighbours): 21 AICc value: 80771.03 
##  Iteration    Log-Likelihood(With bandwidth:  20 )
## =========================
##        0      -1.14e+05 
##        1      -2.556e+06 
##        2      -1.018e+06 
##        3      -3.937e+05 
##        4      -1.649e+05 
##        5      -8.232e+04 
##        6      -5.282e+04 
##        7      -4.313e+04 
##        8      -4.057e+04 
##        9      -4.001e+04 
##       10      -3.989e+04 
##       11      -3.988e+04 
##       12      -3.988e+04 
##       13      -3.988e+04 
## Adaptive bandwidth (number of nearest neighbours): 20 AICc value: 79844.17 
##  Iteration    Log-Likelihood(With bandwidth:  18 )
## =========================
##        0      -1.087e+05 
##        1      -2.794e+06 
##        2      -1.121e+06 
##        3      -4.336e+05 
##        4      -1.793e+05 
##        5      -8.607e+04 
##        6      -5.274e+04 
##        7      -4.193e+04 
##        8      -3.908e+04 
##        9      -3.846e+04 
##       10      -3.832e+04 
##       11      -3.831e+04 
##       12      -3.831e+04 
##       13      -3.831e+04 
## Adaptive bandwidth (number of nearest neighbours): 18 AICc value: 76710.47 
##  Iteration    Log-Likelihood(With bandwidth:  18 )
## =========================
##        0      -1.087e+05 
##        1      -2.794e+06 
##        2      -1.121e+06 
##        3      -4.336e+05 
##        4      -1.793e+05 
##        5      -8.607e+04 
##        6      -5.274e+04 
##        7      -4.193e+04 
##        8      -3.908e+04 
##        9      -3.846e+04 
##       10      -3.832e+04 
##       11      -3.831e+04 
##       12      -3.831e+04 
##       13      -3.831e+04 
## Adaptive bandwidth (number of nearest neighbours): 18 AICc value: 76710.47
gwr <- ggwr.basic(NOR_TXCURR ~ POP_DEN + PHF_COV + AHF_COV + POP_POV + HEALTH_DEP + EDUC_DEP + WORK_DEP + SEC_SHOCK, data = ssdataspdf, bw=bw, family ="poisson", kernel = "exponential", adaptive = FALSE, dMat=DM)
##  Iteration    Log-Likelihood
## =========================
##        0      -3.286e+05 
##        1      -2.356e+11 
##        2      -8.667e+10 
##        3      -3.189e+10 
##        4      -1.173e+10 
##        5      -4.315e+09 
##        6      -1.588e+09 
##        7      -5.844e+08 
##        8      -2.152e+08 
##        9      -7.932e+07 
##       10      -2.931e+07 
##       11      -1.088e+07 
##       12      -4.071e+06 
##       13      -1.548e+06 
##       14      -6.07e+05 
##       15      -2.545e+05 
##       16      -1.242e+05 
##       17      -7.989e+04 
##       18      -6.872e+04 
##       19      -6.726e+04
gwr
##    ***********************************************************************
##    *                       Package   GWmodel                             *
##    ***********************************************************************
##    Program starts at: 2025-01-11 18:20:39 
##    Call:
##    ggwr.basic(formula = NOR_TXCURR ~ POP_DEN + PHF_COV + AHF_COV + 
##     POP_POV + HEALTH_DEP + EDUC_DEP + WORK_DEP + SEC_SHOCK, data = ssdataspdf, 
##     bw = bw, family = "poisson", kernel = "exponential", adaptive = FALSE, 
##     dMat = DM)
## 
##    Dependent (y) variable:  NOR_TXCURR
##    Independent variables:  POP_DEN PHF_COV AHF_COV POP_POV HEALTH_DEP EDUC_DEP WORK_DEP SEC_SHOCK
##    Number of data points: 123
##    Used family: poisson
##    ***********************************************************************
##    *              Results of Generalized linear Regression               *
##    ***********************************************************************
## 
## Call:
## NULL
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
##  -1.000   -0.387   -0.017    0.615  279.945  
## 
## Coefficients:
##                  Estimate     Std. Error z value            Pr(>|z|)    
## Intercept   8.22779520411  0.00628295755 1309.54 <0.0000000000000002 ***
## POP_DEN     0.00006017078  0.00000219293   27.44 <0.0000000000000002 ***
## PHF_COV     0.01568763771  0.00016649451   94.22 <0.0000000000000002 ***
## AHF_COV    -0.00618418686  0.00003826096 -161.63 <0.0000000000000002 ***
## POP_POV    -0.00000317094  0.00000005951  -53.28 <0.0000000000000002 ***
## HEALTH_DEP  0.00000701185  0.00000018534   37.83 <0.0000000000000002 ***
## EDUC_DEP   -0.00003776003  0.00000031498 -119.88 <0.0000000000000002 ***
## WORK_DEP    0.00002833128  0.00000026063  108.70 <0.0000000000000002 ***
## SEC_SHOCK  -0.00006747452  0.00000042839 -157.51 <0.0000000000000002 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 304084  on 122  degrees of freedom
## Residual deviance: 136437  on 114  degrees of freedom
## AIC: 136455
## 
## Number of Fisher Scoring iterations: 5
## 
## 
##  AICc:  136456.5
##  Pseudo R-square value:  0.5513183
##    ***********************************************************************
##    *          Results of Geographically Weighted Regression              *
##    ***********************************************************************
## 
##    *********************Model calibration information*********************
##    Kernel function: exponential 
##    Fixed bandwidth: 18 
##    Regression points: the same locations as observations are used.
##    Distance metric: A distance matrix is specified for this model calibration.
## 
##    ************Summary of Generalized GWR coefficient estimates:**********
##                        Min.       1st Qu.        Median       3rd Qu.    Max.
##    Intercept   8.2111747227  8.2187269022  8.2417798013  8.2560958771  8.2705
##    POP_DEN     0.0000463682  0.0000481389  0.0000516193  0.0000573589  0.0001
##    PHF_COV     0.0149900554  0.0152600207  0.0155711444  0.0157040793  0.0158
##    AHF_COV    -0.0062366651 -0.0061907384 -0.0061425612 -0.0060675648 -0.0060
##    POP_POV    -0.0000033274 -0.0000032653 -0.0000032028 -0.0000030079  0.0000
##    HEALTH_DEP  0.0000068131  0.0000068471  0.0000069027  0.0000069570  0.0000
##    EDUC_DEP   -0.0000372629 -0.0000370117 -0.0000365966 -0.0000362791  0.0000
##    WORK_DEP    0.0000271300  0.0000274056  0.0000282319  0.0000285182  0.0000
##    SEC_SHOCK  -0.0000679106 -0.0000675655 -0.0000673194 -0.0000668842 -0.0001
##    ************************Diagnostic information*************************
##    Number of data points: 123 
##    GW Deviance: 133534.9 
##    AIC : 133554 
##    AICc : 133555.7 
##    Pseudo R-square value:  0.5608618 
## 
##    ***********************************************************************
##    Program stops at: 2025-01-11 18:20:39
bw1 <- bw.ggwr(NOR_TXCURR ~ POP_DEN + PHF_COV + AHF_COV + POP_POV + EDUC_DEP + WORK_DEP + SEC_SHOCK, data = ssdataspdf, family = "poisson", approach = "AICc", kernel = "exponential", adaptive = TRUE, dMat = DM )
##  Iteration    Log-Likelihood(With bandwidth:  83 )
## =========================
##        0      -2.282e+05 
##        1      -1.359e+10 
##        2       -5e+09 
##        3      -1.839e+09 
##        4      -6.768e+08 
##        5      -2.491e+08 
##        6      -9.173e+07 
##        7      -3.382e+07 
##        8      -1.251e+07 
##        9      -4.654e+06 
##       10      -1.753e+06 
##       11      -6.765e+05 
##       12      -2.768e+05 
##       13      -1.292e+05 
##       14      -7.76e+04 
##       15      -6.242e+04 
##       16      -5.903e+04 
##       17      -5.818e+04 
##       18      -5.805e+04 
##       19      -5.805e+04 
## Adaptive bandwidth (number of nearest neighbours): 83 AICc value: 116124.9 
##  Iteration    Log-Likelihood(With bandwidth:  59 )
## =========================
##        0      -1.936e+05 
##        1      -1.654e+08 
##        2      -6.095e+07 
##        3      -2.242e+07 
##        4      -8.265e+06 
##        5      -3.06e+06 
##        6      -1.151e+06 
##        7      -4.467e+05 
##        8      -1.879e+05 
##        9      -9.569e+04 
##       10      -6.538e+04 
##       11      -5.727e+04 
##       12      -5.532e+04 
##       13      -5.485e+04 
##       14      -5.482e+04 
##       15      -5.482e+04 
##       16      -5.482e+04 
## Adaptive bandwidth (number of nearest neighbours): 59 AICc value: 109676.9 
##  Iteration    Log-Likelihood(With bandwidth:  43 )
## =========================
##        0      -1.601e+05 
##        1      -6.701e+06 
##        2      -2.417e+06 
##        3      -8.752e+05 
##        4      -3.323e+05 
##        5      -1.463e+05 
##        6      -8.156e+04 
##        7      -5.892e+04 
##        8      -5.221e+04 
##        9      -5.089e+04 
##       10      -5.066e+04 
##       11      -5.063e+04 
##       12      -5.063e+04 
##       13      -5.063e+04 
## Adaptive bandwidth (number of nearest neighbours): 43 AICc value: 101308.6 
##  Iteration    Log-Likelihood(With bandwidth:  34 )
## =========================
##        0      -1.4e+05 
##        1      -3.842e+06 
##        2      -1.412e+06 
##        3      -5.254e+05 
##        4      -2.153e+05 
##        5      -1.067e+05 
##        6      -6.705e+04 
##        7      -5.288e+04 
##        8      -4.858e+04 
##        9      -4.767e+04 
##       10      -4.75e+04 
##       11      -4.748e+04 
##       12      -4.748e+04 
##       13      -4.748e+04 
## Adaptive bandwidth (number of nearest neighbours): 34 AICc value: 95014.65 
##  Iteration    Log-Likelihood(With bandwidth:  27 )
## =========================
##        0      -1.277e+05 
##        1      -3.227e+06 
##        2      -1.219e+06 
##        3      -4.624e+05 
##        4      -1.919e+05 
##        5      -9.574e+04 
##        6      -6.158e+04 
##        7      -4.946e+04 
##        8      -4.579e+04 
##        9      -4.5e+04 
##       10      -4.484e+04 
##       11      -4.482e+04 
##       12      -4.482e+04 
##       13      -4.482e+04 
## Adaptive bandwidth (number of nearest neighbours): 27 AICc value: 89705.99 
##  Iteration    Log-Likelihood(With bandwidth:  24 )
## =========================
##        0      -1.218e+05 
##        1      -2.85e+06 
##        2      -1.109e+06 
##        3      -4.265e+05 
##        4      -1.778e+05 
##        5      -8.879e+04 
##        6      -5.801e+04 
##        7      -4.737e+04 
##        8      -4.421e+04 
##        9      -4.351e+04 
##       10      -4.336e+04 
##       11      -4.335e+04 
##       12      -4.335e+04 
##       13      -4.335e+04 
## Adaptive bandwidth (number of nearest neighbours): 24 AICc value: 86772.05 
##  Iteration    Log-Likelihood(With bandwidth:  21 )
## =========================
##        0      -1.153e+05 
##        1      -2.433e+06 
##        2      -9.528e+05 
##        3      -3.691e+05 
##        4      -1.558e+05 
##        5      -8.03e+04 
##        6      -5.435e+04 
##        7      -4.538e+04 
##        8      -4.27e+04 
##        9      -4.209e+04 
##       10      -4.196e+04 
##       11      -4.195e+04 
##       12      -4.195e+04 
##       13      -4.195e+04 
## Adaptive bandwidth (number of nearest neighbours): 21 AICc value: 83978.72 
##  Iteration    Log-Likelihood(With bandwidth:  20 )
## =========================
##        0      -1.136e+05 
##        1      -2.228e+06 
##        2      -8.768e+05 
##        3      -3.406e+05 
##        4      -1.445e+05 
##        5      -7.575e+04 
##        6      -5.238e+04 
##        7      -4.443e+04 
##        8      -4.212e+04 
##        9      -4.16e+04 
##       10      -4.149e+04 
##       11      -4.149e+04 
##       12      -4.149e+04 
##       13      -4.149e+04 
## Adaptive bandwidth (number of nearest neighbours): 20 AICc value: 83049.1 
##  Iteration    Log-Likelihood(With bandwidth:  18 )
## =========================
##        0      -1.084e+05 
##        1      -1.998e+06 
##        2      -7.998e+05 
##        3      -3.122e+05 
##        4      -1.321e+05 
##        5      -6.932e+04 
##        6      -4.869e+04 
##        7      -4.22e+04 
##        8      -4.042e+04 
##        9      -4.001e+04 
##       10      -3.994e+04 
##       11      -3.993e+04 
##       12      -3.993e+04 
##       13      -3.993e+04 
## Adaptive bandwidth (number of nearest neighbours): 18 AICc value: 79945.88 
##  Iteration    Log-Likelihood(With bandwidth:  18 )
## =========================
##        0      -1.084e+05 
##        1      -1.998e+06 
##        2      -7.998e+05 
##        3      -3.122e+05 
##        4      -1.321e+05 
##        5      -6.932e+04 
##        6      -4.869e+04 
##        7      -4.22e+04 
##        8      -4.042e+04 
##        9      -4.001e+04 
##       10      -3.994e+04 
##       11      -3.993e+04 
##       12      -3.993e+04 
##       13      -3.993e+04 
## Adaptive bandwidth (number of nearest neighbours): 18 AICc value: 79945.88
gwr1 <- ggwr.basic(NOR_TXCURR ~ POP_DEN + PHF_COV + AHF_COV + POP_POV + HEALTH_DEP + EDUC_DEP + WORK_DEP + SEC_SHOCK, data = ssdataspdf, bw=bw, family ="poisson", kernel = "exponential", adaptive =TRUE, dMat=DM)
##  Iteration    Log-Likelihood
## =========================
##        0      -1.087e+05 
##        1      -2.794e+06 
##        2      -1.121e+06 
##        3      -4.336e+05 
##        4      -1.793e+05 
##        5      -8.607e+04 
##        6      -5.274e+04 
##        7      -4.193e+04 
##        8      -3.908e+04 
##        9      -3.846e+04 
##       10      -3.832e+04 
##       11      -3.831e+04 
##       12      -3.831e+04 
##       13      -3.831e+04
gwr1
##    ***********************************************************************
##    *                       Package   GWmodel                             *
##    ***********************************************************************
##    Program starts at: 2025-01-11 18:20:40 
##    Call:
##    ggwr.basic(formula = NOR_TXCURR ~ POP_DEN + PHF_COV + AHF_COV + 
##     POP_POV + HEALTH_DEP + EDUC_DEP + WORK_DEP + SEC_SHOCK, data = ssdataspdf, 
##     bw = bw, family = "poisson", kernel = "exponential", adaptive = TRUE, 
##     dMat = DM)
## 
##    Dependent (y) variable:  NOR_TXCURR
##    Independent variables:  POP_DEN PHF_COV AHF_COV POP_POV HEALTH_DEP EDUC_DEP WORK_DEP SEC_SHOCK
##    Number of data points: 123
##    Used family: poisson
##    ***********************************************************************
##    *              Results of Generalized linear Regression               *
##    ***********************************************************************
## 
## Call:
## NULL
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
##  -1.000   -0.387   -0.017    0.615  279.945  
## 
## Coefficients:
##                  Estimate     Std. Error z value            Pr(>|z|)    
## Intercept   8.22779520411  0.00628295755 1309.54 <0.0000000000000002 ***
## POP_DEN     0.00006017078  0.00000219293   27.44 <0.0000000000000002 ***
## PHF_COV     0.01568763771  0.00016649451   94.22 <0.0000000000000002 ***
## AHF_COV    -0.00618418686  0.00003826096 -161.63 <0.0000000000000002 ***
## POP_POV    -0.00000317094  0.00000005951  -53.28 <0.0000000000000002 ***
## HEALTH_DEP  0.00000701185  0.00000018534   37.83 <0.0000000000000002 ***
## EDUC_DEP   -0.00003776003  0.00000031498 -119.88 <0.0000000000000002 ***
## WORK_DEP    0.00002833128  0.00000026063  108.70 <0.0000000000000002 ***
## SEC_SHOCK  -0.00006747452  0.00000042839 -157.51 <0.0000000000000002 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 304084  on 122  degrees of freedom
## Residual deviance: 136437  on 114  degrees of freedom
## AIC: 136455
## 
## Number of Fisher Scoring iterations: 5
## 
## 
##  AICc:  136456.5
##  Pseudo R-square value:  0.5513183
##    ***********************************************************************
##    *          Results of Geographically Weighted Regression              *
##    ***********************************************************************
## 
##    *********************Model calibration information*********************
##    Kernel function: exponential 
##    Adaptive bandwidth: 18 (number of nearest neighbours)
##    Regression points: the same locations as observations are used.
##    Distance metric: A distance matrix is specified for this model calibration.
## 
##    ************Summary of Generalized GWR coefficient estimates:**********
##                        Min.       1st Qu.        Median       3rd Qu.    Max.
##    Intercept   7.5694535226  7.8182781956  8.1580282754  8.5895840674  9.2153
##    POP_DEN    -0.0001761390 -0.0000718419  0.0000038045  0.0001266410  0.0003
##    PHF_COV    -0.0293228488  0.0061415428  0.0137416914  0.0213268175  0.0291
##    AHF_COV    -0.0116756893 -0.0083821354 -0.0059160307 -0.0046285705 -0.0032
##    POP_POV    -0.0000070754 -0.0000048468 -0.0000039438 -0.0000007892  0.0000
##    HEALTH_DEP -0.0000434281 -0.0000064849  0.0000054519  0.0000080731  0.0000
##    EDUC_DEP   -0.0000423468 -0.0000273992 -0.0000230528 -0.0000082143  0.0001
##    WORK_DEP    0.0000021100  0.0000118604  0.0000317444  0.0000430481  0.0001
##    SEC_SHOCK  -0.0000754347 -0.0000660665 -0.0000502261 -0.0000392972  0.0000
##    ************************Diagnostic information*************************
##    Number of data points: 123 
##    GW Deviance: 75629.48 
##    AIC : 75692.78 
##    AICc : 75715.66 
##    Pseudo R-square value:  0.7512876 
## 
##    ***********************************************************************
##    Program stops at: 2025-01-11 18:20:40
gwr$bandwidth
## NULL
gwr$SDF
## class       : SpatialPolygonsDataFrame 
## features    : 123 
## extent      : 4.978159, 9.467367, 4.270418, 7.573423  (xmin, xmax, ymin, ymax)
## crs         : +proj=longlat +datum=WGS84 +no_defs 
## variables   : 30
## names       :       Intercept,               POP_DEN,            PHF_COV,              AHF_COV,                 POP_POV,             HEALTH_DEP,               EDUC_DEP,              WORK_DEP,              SEC_SHOCK,     y,              yhat,          residual,       Intercept_SE,            POP_DEN_SE,           PHF_COV_SE, ... 
## min values  : 8.2111747226557, 0.0000463682087169116, 0.0149900553967944, -0.00623666506062062, -0.00000332739638422906,  0.0000068130586121287,  -0.000037262860410575, 0.0000271300062146162, -0.0000679106363728112,     0, 0.498680455500644, -4845.78772155221, 0.0443679751755873, 0.0000108111713605242, 0.000313542139110611, ... 
## max values  : 8.2704913487846, 0.0000615754326115069, 0.0158318110014362, -0.00604208315301457, -0.00000289996024892797, 0.00000703381128054837, -0.0000358893557158914, 0.0000288739846652074, -0.0000665217652478589, 13299,  7781.74896182731,  6429.20407371782, 0.0467326873032881, 0.0000112903996107685, 0.000342760162663926, ...
gc <- round(cor(as.data.frame(gwr$SDF[,2:11]), use ="complete.obs"),2)
corrplot.mixed(gc, tl.cex = 0.8, tl.col="black")

pairs(as(gwr$SDF, "data.frame")[,2:11], pch=".")

gwr$GW.diagnostic
## $gw.deviance
##        1 
## 133534.9 
## 
## $AICc
##        1 
## 133555.7 
## 
## $AIC
##      1 
## 133554 
## 
## $pseudo.R2
##         1 
## 0.5608618 
## 
## $edf
## [1] 14866.85

Create spatial data frame

ssdataspdf@data$y<-gwr$SDF$y
ssdataspdf@data$yhat<-gwr$SDF$yhat
ssdataspdf@data$residual<-gwr$SDF$residual
rsd=sd(ssdataspdf@data$residual)
ssdataspdf@data$stdRes<-(ssdataspdf@data$residual)/sd(ssdataspdf@data$residual)
ssdataspdf@data$LLN=ssdataspdf@data$yhat-1.645*rsd
ssdataspdf@data$ULN=ssdataspdf@data$yhat+1.645*rsd



# Intercept
ssdataspdf@data$Intercept<-gwr$SDF$Intercept
ssdataspdf@data$est_POP_DEN<-gwr$SDF$POP_DEN
ssdataspdf@data$est_PHF_COV<-gwr$SDF$PHF_COV
ssdataspdf@data$est_AHF_COV<-gwr$SDF$AHF_COV
ssdataspdf@data$est_POP_POV<-gwr$SDF$POP_POV
ssdataspdf@data$est_HEALTH_DEP<-gwr$SDF$HEALTH_DEP
ssdataspdf@data$est_EDUC_DEP<-gwr$SDF$EDUC_DEP
ssdataspdf@data$est_WORK_DEP<-gwr$SDF$WORK_DEP
ssdataspdf@data$est_SEC_SHOCK<-gwr$SDF$SEC_SHOCK

# T-values
ssdataspdf@data$t_Intercept<-gwr$SDF$Intercept_TV
ssdataspdf@data$t_POP_DEN<-gwr$SDF$POP_DEN_TV
ssdataspdf@data$t_PHF_COV<-gwr$SDF$PHF_COV_TV
ssdataspdf@data$t_AHF_COV<-gwr$SDF$AHF_COV_TV
ssdataspdf@data$t_POP_POV<-gwr$SDF$POP_POV_TV
ssdataspdf@data$t_HEALTH_DEP<-gwr$SDF$HEALTH_DEP_TV
ssdataspdf@data$t_EDUC_DEP<-gwr$SDF$EDUC_DEP_TV
ssdataspdf@data$t_WORK_DEP<-gwr$SDF$WORK_DEP_TV
ssdataspdf@data$t_SEC_SHOCK<-gwr$SDF$SEC_SHOCK_TV


# Calculate psudo-t values

ssdataspdf@data$p_POP_DEN<-2*pt(-abs(gwr$SDF$POP_DEN_TV),df=3103)
ssdataspdf@data$p_PHF_COV<-2*pt(-abs(gwr$SDF$PHF_COV_TV),df=3103)
ssdataspdf@data$p_AHF_COV<-2*pt(-abs(gwr$SDF$AHF_COV_TV),df=3103)
ssdataspdf@data$p_POP_POV<-2*pt(-abs(gwr$SDF$POP_POV_TV),df=3103)
ssdataspdf@data$p_HEALTH_DEP<-2*pt(-abs(gwr$SDF$HEALTH_DEP_TV),df=3103)
ssdataspdf@data$p_EDUC_DEP<-2*pt(-abs(gwr$SDF$EDUC_DEP_TV),df=3103)
ssdataspdf@data$p_WORK_DEP<-2*pt(-abs(gwr$SDF$WORK_DEP_TV),df=3103)
ssdataspdf@data$p_SEC_SHOCK<-2*pt(-abs(gwr$SDF$SEC_SHOCK_TV),df=3103)

ssdataspdf$sig_POP_DEN <-ifelse(ssdataspdf@data$est_POP_DEN > 0 &
                               ssdataspdf@data$p_POP_DEN <= 0.05 , 'Yes', "No")
ssdataspdf$sig_PHF_COV <-ifelse(ssdataspdf@data$est_PHF_COV > 0 &
                              ssdataspdf@data$p_PHF_COV <= 0.05 , 'Yes', "No")
ssdataspdf$sig_AHF_COV <-ifelse(ssdataspdf@data$est_AHF_COV > 0 &
                               ssdataspdf@data$p_AHF_COV <= 0.05 , 'Yes', "No")
ssdataspdf$sig_POP_POV <-ifelse(ssdataspdf@data$est_POP_POV > 0 &
                              ssdataspdf@data$p_POP_POV <= 0.05 , 'Yes', "No")
ssdataspdf$sig_HEALTH_DEP <-ifelse(ssdataspdf@data$est_HEALTH_DEP > 0 &
                              ssdataspdf@data$p_HEALTH_DEP <= 0.05 , 'Yes', "No")
ssdataspdf$sig_EDUC_DEP <-ifelse(ssdataspdf@data$est_EDUC_DEP > 0 &
                                     ssdataspdf@data$p_EDUC_DEP <= 0.05 , 'Yes', "No")
ssdataspdf$sig_WORK_DEP <-ifelse(ssdataspdf@data$est_WORK_DEP > 0 &
                                     ssdataspdf@data$p_WORK_DEP <= 0.05 , 'Yes', "No")
ssdataspdf$sig_SEC_SHOCK <-ifelse(ssdataspdf@data$est_SEC_SHOCK > 0 &
                                     ssdataspdf@data$p_SEC_SHOCK <= 0.05 , 'Yes', "No")

#Local estimates

popden <- tm_shape(ssdataspdf) + tm_polygons(col= "est_POP_DEN", palette = "GnBu", border.col = "grey80", border.alpha = 0.5, title = "", style = "cont", legend.is.portrait = FALSE) + tm_layout(legend.outside.position = "bottom", legend.outside.size = 0.2, legend.text.size = 5, legend.outside = TRUE, main.title = "Population Density", main.title.size = 0.75) 
phfcov <- tm_shape(ssdataspdf) + tm_polygons(col= "est_PHF_COV", palette = "GnBu", border.col = "grey80", border.alpha = 0.5, title = "", style = "cont", legend.is.portrait = FALSE) + tm_layout(legend.outside.position = "bottom", legend.outside.size = 0.2, legend.text.size = 5, legend.outside = TRUE, main.title = "PHF Coverage", main.title.size = 0.75) 
ahfcov <- tm_shape(ssdataspdf) + tm_polygons(col= "est_AHF_COV", palette = "GnBu", border.col = "grey80", border.alpha = 0.5, title = "", style = "cont", legend.is.portrait = FALSE) + tm_layout(legend.outside.position = "bottom", legend.outside.size = 0.2, legend.text.size = 5, legend.outside = TRUE, main.title = "ART HF Coverage", main.title.size = 0.75) 
ppov <- tm_shape(ssdataspdf) + tm_polygons(col= "est_POP_POV", palette = "GnBu", border.col = "grey80", border.alpha = 0.5, title = "", style = "cont", legend.is.portrait = FALSE) + tm_layout(legend.outside.position = "bottom", legend.outside.size = 0.2, legend.text.size = 5, legend.outside = TRUE, main.title = "Population Poverty", main.title.size = 0.75) 
hdev <- tm_shape(ssdataspdf) + tm_polygons(col= "est_HEALTH_DEP", palette = "GnBu", border.col = "grey80", border.alpha = 0.5, title = "", style = "cont", legend.is.portrait = FALSE) + tm_layout(legend.outside.position = "bottom", legend.outside.size = 0.2, legend.text.size = 5, legend.outside = TRUE, main.title = "Health Deprivation", main.title.size = 0.75) 
edev <- tm_shape(ssdataspdf) + tm_polygons(col= "est_EDUC_DEP", palette = "GnBu", border.col = "grey80", border.alpha = 0.5, title = "", style = "cont", legend.is.portrait = FALSE) + tm_layout(legend.outside.position = "bottom", legend.outside.size = 0.2, legend.text.size = 5, legend.outside = TRUE, main.title = "Education Deprivation", main.title.size = 0.75) 
wdev <- tm_shape(ssdataspdf) + tm_polygons(col= "est_WORK_DEP", palette = "GnBu", border.col = "grey80", border.alpha = 0.5, title = "", style = "cont", legend.is.portrait = FALSE) + tm_layout(legend.outside.position = "bottom", legend.outside.size = 0.2, legend.text.size = 5, legend.outside = TRUE, main.title = "Work Deprivation", main.title.size = 0.75) 
sshk <- tm_shape(ssdataspdf) + tm_polygons(col= "est_SEC_SHOCK", palette = "GnBu", border.col = "grey80", border.alpha = 0.5, title = "", style = "cont", legend.is.portrait = FALSE) + tm_layout(legend.outside.position = "bottom", legend.outside.size = 0.2, legend.text.size = 5, legend.outside = TRUE, main.title = "Security Shock", main.title.size = 0.75) 
tmap_arrange(popden, phfcov, ahfcov, ppov, hdev, edev, wdev, sshk, ncol = 4)
## Legend labels were too wide. Therefore, legend.text.size has been set to 0.64. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.
## Legend labels were too wide. Therefore, legend.text.size has been set to 0.51. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.
## Legend labels were too wide. Therefore, legend.text.size has been set to 0.5. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.
## Legend labels were too wide. Therefore, legend.text.size has been set to 0.32. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.
## Legend labels were too wide. Therefore, legend.text.size has been set to 0.38. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.
## Legend labels were too wide. Therefore, legend.text.size has been set to 0.22. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.
## Legend labels were too wide. Therefore, legend.text.size has been set to 0.56. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.
## Legend labels were too wide. Therefore, legend.text.size has been set to 0.22. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.

tpopden <- tm_shape(ssdataspdf) + tm_polygons(col= "t_POP_DEN", palette = "GnBu", border.col = "grey80", border.alpha = 0.5, title = "", style = "cont", legend.is.portrait = FALSE) + tm_layout(legend.outside.position = "bottom", legend.outside.size = 0.2, legend.text.size = 5, legend.outside = TRUE, main.title = "Population Density", main.title.size = 0.75) 
tphfcov <- tm_shape(ssdataspdf) + tm_polygons(col= "t_PHF_COV", palette = "GnBu", border.col = "grey80", border.alpha = 0.5, title = "", style = "cont", legend.is.portrait = FALSE) + tm_layout(legend.outside.position = "bottom", legend.outside.size = 0.2, legend.text.size = 5, legend.outside = TRUE, main.title = "PHF Coverage", main.title.size = 0.75) 
tahfcov <- tm_shape(ssdataspdf) + tm_polygons(col= "t_AHF_COV", palette = "GnBu", border.col = "grey80", border.alpha = 0.5, title = "", style = "cont", legend.is.portrait = FALSE) + tm_layout(legend.outside.position = "bottom", legend.outside.size = 0.2, legend.text.size = 5, legend.outside = TRUE, main.title = "ART HF Coverage", main.title.size = 0.75) 
tppov <- tm_shape(ssdataspdf) + tm_polygons(col= "t_POP_POV", palette = "GnBu", border.col = "grey80", border.alpha = 0.5, title = "", style = "cont", legend.is.portrait = FALSE) + tm_layout(legend.outside.position = "bottom", legend.outside.size = 0.2, legend.text.size = 5, legend.outside = TRUE, main.title = "Population Poverty", main.title.size = 0.75) 
thdev <- tm_shape(ssdataspdf) + tm_polygons(col= "t_HEALTH_DEP", palette = "GnBu", border.col = "grey80", border.alpha = 0.5, title = "", style = "cont", legend.is.portrait = FALSE) + tm_layout(legend.outside.position = "bottom", legend.outside.size = 0.2, legend.text.size = 5, legend.outside = TRUE, main.title = "Health Deprivation", main.title.size = 0.75) 
tedev <- tm_shape(ssdataspdf) + tm_polygons(col= "t_EDUC_DEP", palette = "GnBu", border.col = "grey80", border.alpha = 0.5, title = "", style = "cont", legend.is.portrait = FALSE) + tm_layout(legend.outside.position = "bottom", legend.outside.size = 0.2, legend.text.size = 5, legend.outside = TRUE, main.title = "Education Deprivation", main.title.size = 0.75) 
twdev <- tm_shape(ssdataspdf) + tm_polygons(col= "t_WORK_DEP", palette = "GnBu", border.col = "grey80", border.alpha = 0.5, title = "", style = "cont", legend.is.portrait = FALSE) + tm_layout(legend.outside.position = "bottom", legend.outside.size = 0.2, legend.text.size = 5, legend.outside = TRUE, main.title = "Work Deprivation", main.title.size = 0.75) 
tsshk <- tm_shape(ssdataspdf) + tm_polygons(col= "t_SEC_SHOCK", palette = "GnBu", border.col = "grey80", border.alpha = 0.5, title = "", style = "cont", legend.is.portrait = FALSE) + tm_layout(legend.outside.position = "bottom", legend.outside.size = 0.2, legend.text.size = 5, legend.outside = TRUE, main.title = "Security Shock", main.title.size = 0.75) 
tmap_arrange(tpopden, tphfcov, tahfcov, tppov, thdev, tedev, twdev, tsshk, ncol = 4)
## Legend labels were too wide. Therefore, legend.text.size has been set to 0.64. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.
## Legend labels were too wide. Therefore, legend.text.size has been set to 0.55. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.

spopden <- tm_shape(ssdataspdf) + tm_polygons(col= "sig_POP_DEN", palette = "GnBu", border.col = "grey80", border.alpha = 0.5, title = "", style = "cont", legend.is.portrait = FALSE) + tm_layout(legend.outside.position = "bottom", legend.outside.size = 0.2, legend.text.size = 5, legend.outside = TRUE, main.title = "Population Density", main.title.size = 0.75) 
sphfcov <- tm_shape(ssdataspdf) + tm_polygons(col= "sig_PHF_COV", palette = "GnBu", border.col = "grey80", border.alpha = 0.5, title = "", style = "cont", legend.is.portrait = FALSE) + tm_layout(legend.outside.position = "bottom", legend.outside.size = 0.2, legend.text.size = 5, legend.outside = TRUE, main.title = "PHF Coverage", main.title.size = 0.75) 
sahfcov <- tm_shape(ssdataspdf) + tm_polygons(col= "sig_AHF_COV", palette = "GnBu", border.col = "grey80", border.alpha = 0.5, title = "", style = "cont", legend.is.portrait = FALSE) + tm_layout(legend.outside.position = "bottom", legend.outside.size = 0.2, legend.text.size = 5, legend.outside = TRUE, main.title = "ART HF Coverage", main.title.size = 0.75) 
sppov <- tm_shape(ssdataspdf) + tm_polygons(col= "sig_POP_POV", palette = "GnBu", border.col = "grey80", border.alpha = 0.5, title = "", style = "cont", legend.is.portrait = FALSE) + tm_layout(legend.outside.position = "bottom", legend.outside.size = 0.2, legend.text.size = 5, legend.outside = TRUE, main.title = "Population Poverty", main.title.size = 0.75) 
shdev <- tm_shape(ssdataspdf) + tm_polygons(col= "sig_HEALTH_DEP", palette = "GnBu", border.col = "grey80", border.alpha = 0.5, title = "", style = "cont", legend.is.portrait = FALSE) + tm_layout(legend.outside.position = "bottom", legend.outside.size = 0.2, legend.text.size = 5, legend.outside = TRUE, main.title = "Health Deprivation", main.title.size = 0.75) 
sedev <- tm_shape(ssdataspdf) + tm_polygons(col= "sig_EDUC_DEP", palette = "GnBu", border.col = "grey80", border.alpha = 0.5, title = "", style = "cont", legend.is.portrait = FALSE) + tm_layout(legend.outside.position = "bottom", legend.outside.size = 0.2, legend.text.size = 5, legend.outside = TRUE, main.title = "Education Deprivation", main.title.size = 0.75) 
swdev <- tm_shape(ssdataspdf) + tm_polygons(col= "sig_WORK_DEP", palette = "GnBu", border.col = "grey80", border.alpha = 0.5, title = "", style = "cont", legend.is.portrait = FALSE) + tm_layout(legend.outside.position = "bottom", legend.outside.size = 0.2, legend.text.size = 0.55, legend.outside = TRUE, main.title = "Work Deprivation", main.title.size = 0.75) 
ssshk <- tm_shape(ssdataspdf) + tm_polygons(col= "sig_SEC_SHOCK", palette = "GnBu", border.col = "grey80", border.alpha = 0.5, title = "", style = "cont", legend.is.portrait = FALSE) + tm_layout(legend.outside.position = "bottom", legend.outside.size = 0.2, legend.text.size = 5, legend.outside = TRUE, main.title = "Security Shock", main.title.size = 0.75) 
tmap_arrange(spopden, sphfcov, sahfcov, sppov, shdev, sedev, swdev, ssshk, ncol = 4)

tm_shape(ssdataspdf) + tm_polygons(col= "stdRes", palette = "GnBu", border.col = "grey80", border.alpha = 0.5, title = "GWRP Std. Residuals", style = "cont")

library(writexl)
## Warning: package 'writexl' was built under R version 4.2.3
write_xlsx(ssdataspdf@data, "C:/Users/xyaze/Desktop/ssdataspdf.xls")
library(spatialreg)
## Warning: package 'spatialreg' was built under R version 4.2.3
## Loading required package: Matrix
## 
## Attaching package: 'spatialreg'
## The following objects are masked from 'package:spdep':
## 
##     get.ClusterOption, get.coresOption, get.mcOption,
##     get.VerboseOption, get.ZeroPolicyOption, set.ClusterOption,
##     set.coresOption, set.mcOption, set.VerboseOption,
##     set.ZeroPolicyOption
ev <- eigenw(sscastlistw)
W <- as(sscastlistw, "CsparseMatrix")
trMatc <- trW(W, type="mult")
txlageig <- lagsarlm(NOR_TXCURR ~ POP_DEN + PHF_COV + AHF_COV + POP_POV + HEALTH_DEP + EDUC_DEP + WORK_DEP + SEC_SHOCK,
data = rschdata, listw=sscastlistw,
method="eigen", quiet=FALSE, control=list(pre_eig=ev, OrdVsign=1))
## 
## Spatial lag model
## Jacobian calculated using neighbourhood matrix eigenvalues
## 
## rho:  -0.5735677     function value:  -1146.661 
## rho:  0.02748166     function value:  -1109.931 
## rho:  0.3989506  function value:  -1095.66 
## rho:  0.628531   function value:  -1094.018 
## rho:  0.5824889  function value:  -1093.742 
## rho:  0.5636673  function value:  -1093.725 
## rho:  0.5674025  function value:  -1093.724 
## rho:  0.5671888  function value:  -1093.724 
## rho:  0.5672039  function value:  -1093.724 
## rho:  0.5672041  function value:  -1093.724 
## rho:  0.5672041  function value:  -1093.724 
## rho:  0.567204   function value:  -1093.724 
## rho:  0.567204   function value:  -1093.724 
## rho:  0.567204   function value:  -1093.724 
## rho:  0.567204   function value:  -1093.724
## Warning in lagsarlm(NOR_TXCURR ~ POP_DEN + PHF_COV + AHF_COV + POP_POV + : inversion of asymptotic covariance matrix failed for tol.solve = 0.000000000000000222044604925031 
##   reciprocal condition number = 2.55991e-18 - using numerical Hessian.
## Hessian: rho:     0.567204   function value:  -1093.724 
## Hessian: rho:     0.5672075  function value:  -1093.724 
## Hessian: rho:     0.567204   function value:  -1093.724 
## Hessian: rho:     0.567204   function value:  -1093.724 
## Hessian: rho:     0.567204   function value:  -1093.724 
## Hessian: rho:     0.567204   function value:  -1093.724 
## Hessian: rho:     0.567204   function value:  -1093.724 
## Hessian: rho:     0.567204   function value:  -1093.724 
## Hessian: rho:     0.567204   function value:  -1093.724 
## Hessian: rho:     0.567204   function value:  -1093.724 
## Hessian: rho:     0.567204   function value:  -1093.724 
## Hessian: rho:     0.5672006  function value:  -1093.724 
## Hessian: rho:     0.567204   function value:  -1093.724 
## Hessian: rho:     0.567204   function value:  -1093.724 
## Hessian: rho:     0.567204   function value:  -1093.724 
## Hessian: rho:     0.567204   function value:  -1093.724 
## Hessian: rho:     0.567204   function value:  -1093.724 
## Hessian: rho:     0.567204   function value:  -1093.724 
## Hessian: rho:     0.567204   function value:  -1093.724 
## Hessian: rho:     0.567204   function value:  -1093.724 
## Hessian: rho:     0.567204   function value:  -1093.724 
## Hessian: rho:     0.5672075  function value:  -1093.724 
## Hessian: rho:     0.5672075  function value:  -1093.724 
## Hessian: rho:     0.5672075  function value:  -1093.724 
## Hessian: rho:     0.5672075  function value:  -1093.724 
## Hessian: rho:     0.5672075  function value:  -1093.724 
## Hessian: rho:     0.5672075  function value:  -1093.724 
## Hessian: rho:     0.5672075  function value:  -1093.724 
## Hessian: rho:     0.5672075  function value:  -1093.724 
## Hessian: rho:     0.5672075  function value:  -1093.724 
## Hessian: rho:     0.567204   function value:  -1093.724 
## Hessian: rho:     0.567204   function value:  -1093.724 
## Hessian: rho:     0.567204   function value:  -1093.724 
## Hessian: rho:     0.567204   function value:  -1093.724 
## Hessian: rho:     0.567204   function value:  -1093.724 
## Hessian: rho:     0.567204   function value:  -1093.724 
## Hessian: rho:     0.567204   function value:  -1093.724 
## Hessian: rho:     0.567204   function value:  -1093.724 
## Hessian: rho:     0.567204   function value:  -1093.724 
## Hessian: rho:     0.567204   function value:  -1093.724 
## Hessian: rho:     0.567204   function value:  -1093.724 
## Hessian: rho:     0.567204   function value:  -1093.724 
## Hessian: rho:     0.567204   function value:  -1093.724 
## Hessian: rho:     0.567204   function value:  -1093.724 
## Hessian: rho:     0.567204   function value:  -1093.724 
## Hessian: rho:     0.567204   function value:  -1093.724 
## Hessian: rho:     0.567204   function value:  -1093.724 
## Hessian: rho:     0.567204   function value:  -1093.724 
## Hessian: rho:     0.567204   function value:  -1093.724 
## Hessian: rho:     0.567204   function value:  -1093.724 
## Hessian: rho:     0.567204   function value:  -1093.724 
## Hessian: rho:     0.567204   function value:  -1093.724 
## Hessian: rho:     0.567204   function value:  -1093.724 
## Hessian: rho:     0.567204   function value:  -1093.724 
## Hessian: rho:     0.567204   function value:  -1093.724 
## Hessian: rho:     0.567204   function value:  -1093.724 
## Hessian: rho:     0.567204   function value:  -1093.724 
## Hessian: rho:     0.567204   function value:  -1093.724 
## Hessian: rho:     0.567204   function value:  -1093.724 
## Hessian: rho:     0.567204   function value:  -1093.724 
## Hessian: rho:     0.567204   function value:  -1093.724 
## Hessian: rho:     0.567204   function value:  -1093.724 
## Hessian: rho:     0.567204   function value:  -1093.724 
## Hessian: rho:     0.567204   function value:  -1093.724 
## Hessian: rho:     0.567204   function value:  -1093.724 
## Hessian: rho:     0.567204   function value:  -1093.724
txx <- summary(txlageig, correlation=TRUE)
coef(txx)
##                   Estimate    Std. Error   z value    Pr(>|z|)
## (Intercept) 1341.926107743 561.250367064  2.390958 0.016804471
## POP_DEN        0.173660156   0.211841653  0.819764 0.412350670
## PHF_COV       42.583424408  16.504786307  2.580065 0.009878167
## AHF_COV       -1.762526576   0.828863099 -2.126439 0.033466748
## POP_POV       -0.007285443   0.004060498 -1.794224 0.072777412
## HEALTH_DEP     0.021234775   0.013142936  1.615680 0.106163586
## EDUC_DEP      -0.064440087   0.021990504 -2.930360 0.003385698
## WORK_DEP       0.014471621   0.014182363  1.020396 0.307540858
## SEC_SHOCK     -0.050739839   0.026472907 -1.916670 0.055279806
txx
## 
## Call:lagsarlm(formula = NOR_TXCURR ~ POP_DEN + PHF_COV + AHF_COV + 
##     POP_POV + HEALTH_DEP + EDUC_DEP + WORK_DEP + SEC_SHOCK, data = rschdata, 
##     listw = sscastlistw, method = "eigen", quiet = FALSE, control = list(pre_eig = ev, 
##         OrdVsign = 1))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -3050.92  -890.83  -172.36   586.74  7774.97 
## 
## Type: lag 
## Coefficients: (numerical Hessian approximate standard errors) 
##                 Estimate   Std. Error z value Pr(>|z|)
## (Intercept) 1341.9261077  561.2503671  2.3910 0.016804
## POP_DEN        0.1736602    0.2118417  0.8198 0.412351
## PHF_COV       42.5834244   16.5047863  2.5801 0.009878
## AHF_COV       -1.7625266    0.8288631 -2.1264 0.033467
## POP_POV       -0.0072854    0.0040605 -1.7942 0.072777
## HEALTH_DEP     0.0212348    0.0131429  1.6157 0.106164
## EDUC_DEP      -0.0644401    0.0219905 -2.9304 0.003386
## WORK_DEP       0.0144716    0.0141824  1.0204 0.307541
## SEC_SHOCK     -0.0507398    0.0264729 -1.9167 0.055280
## 
## Rho: 0.5672, LR test value: 35.26, p-value: 0.0000000028847
## Approximate (numerical Hessian) standard error: 0.081522
##     z-value: 6.9577, p-value: 0.0000000000034583
## Wald statistic: 48.41, p-value: 0.0000000000034583
## 
## Log likelihood: -1093.724 for lag model
## ML residual variance (sigma squared): 2866000, (sigma: 1692.9)
## Number of observations: 123 
## Number of parameters estimated: 11 
## AIC: 2209.4, (AIC for lm: 2242.7)
## 
##  Approximate correlation of coefficients 
##             rho   (Intercept) POP_DEN PHF_COV AHF_COV POP_POV HEALTH_DEP
## (Intercept) -0.36                                                       
## POP_DEN     -0.14 -0.14                                                 
## PHF_COV      0.08 -0.70        0.12                                     
## AHF_COV      0.27 -0.31        0.04    0.03                             
## POP_POV     -0.05  0.07        0.31   -0.06   -0.16                     
## HEALTH_DEP   0.02 -0.25       -0.33    0.14    0.21   -0.70             
## EDUC_DEP     0.07 -0.02        0.33   -0.01   -0.12    0.70   -0.70     
## WORK_DEP    -0.23  0.16       -0.15    0.02   -0.12   -0.17   -0.33     
## SEC_SHOCK    0.25 -0.18        0.04    0.01    0.01   -0.22   -0.16     
##             EDUC_DEP WORK_DEP
## (Intercept)                  
## POP_DEN                      
## PHF_COV                      
## AHF_COV                      
## POP_POV                      
## HEALTH_DEP                   
## EDUC_DEP                     
## WORK_DEP    -0.12            
## SEC_SHOCK   -0.15    -0.05
txerrorlag <- errorsarlm(NOR_TXCURR ~ POP_DEN + PHF_COV + AHF_COV + POP_POV + HEALTH_DEP + EDUC_DEP + WORK_DEP + SEC_SHOCK,
data = rschdata, listw=sscastlistw)
## Warning in errorsarlm(NOR_TXCURR ~ POP_DEN + PHF_COV + AHF_COV + POP_POV + : inversion of asymptotic covariance matrix failed for tol.solve = 0.000000000000000222044604925031 
##   reciprocal condition number = 8.18826e-18 - using numerical Hessian.
txerrorlag
## 
## Call:
## errorsarlm(formula = NOR_TXCURR ~ POP_DEN + PHF_COV + AHF_COV + 
##     POP_POV + HEALTH_DEP + EDUC_DEP + WORK_DEP + SEC_SHOCK, data = rschdata, 
##     listw = sscastlistw)
## Type: error 
## 
## Coefficients:
##         lambda    (Intercept)        POP_DEN        PHF_COV        AHF_COV 
##    0.702144460 2288.142081499    0.086918521   48.297196678   -0.346993647 
##        POP_POV     HEALTH_DEP       EDUC_DEP       WORK_DEP      SEC_SHOCK 
##   -0.008551277    0.008546915   -0.049243440    0.017489395    0.001965211 
## 
## Log likelihood: -1097.267
summary(txerrorlag)
## 
## Call:errorsarlm(formula = NOR_TXCURR ~ POP_DEN + PHF_COV + AHF_COV + 
##     POP_POV + HEALTH_DEP + EDUC_DEP + WORK_DEP + SEC_SHOCK, data = rschdata, 
##     listw = sscastlistw)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -3013.85  -932.99  -161.60   482.54  7734.31 
## 
## Type: error 
## Coefficients: (asymptotic standard errors) 
##                 Estimate   Std. Error z value  Pr(>|z|)
## (Intercept) 2288.1420815  689.5807901  3.3182 0.0009061
## POP_DEN        0.0869185    0.1972520  0.4406 0.6594685
## PHF_COV       48.2971967   15.5461510  3.1067 0.0018919
## AHF_COV       -0.3469936    0.8161494 -0.4252 0.6707204
## POP_POV       -0.0085513    0.0038551 -2.2182 0.0265426
## HEALTH_DEP     0.0085469    0.0124485  0.6866 0.4923476
## EDUC_DEP      -0.0492434    0.0196173 -2.5102 0.0120660
## WORK_DEP       0.0174894    0.0175883  0.9944 0.3200407
## SEC_SHOCK      0.0019652    0.0291219  0.0675 0.9461979
## 
## Lambda: 0.70214, LR test value: 28.176, p-value: 0.00000011079
## Approximate (numerical Hessian) standard error: 0.076702
##     z-value: 9.1542, p-value: < 0.000000000000000222
## Wald statistic: 83.8, p-value: < 0.000000000000000222
## 
## Log likelihood: -1097.267 for error model
## ML residual variance (sigma squared): 2877300, (sigma: 1696.3)
## Number of observations: 123 
## Number of parameters estimated: 11 
## AIC: NA (not available for weighted model), (AIC for lm: 2242.7)
ssdata$pos <- predict(poisson_output)
ssdata$nbr <- predict(gnb)
ssdata$txeig <- predict(txlageig)
## This method assumes the response is known - see manual page
ssdata$txerr <- predict(txerrorlag)
## This method assumes the response is known - see manual page
tm_shape(ssdata) + tm_polygons(col = "nbr", palette = "GnBu", title = "", border.col = "grey80", border.alpha = 0.5)

tm_shape(ssdata) + tm_polygons(col = "pos", palette = "GnBu", title = "", border.col = "grey80", border.alpha = 0.5)

tm_shape(ssdata) + tm_polygons(col = "txeig", palette = "GnBu", title = "txeig", border.col = "grey80", border.alpha = 0.5)

tm_shape(ssdata) + tm_polygons(col = "txerr", palette = "GnBu", title = "txerr", border.col = "grey80", border.alpha = 0.5)

ssdata$posres <- residuals(poisson_output)
ssdata$nbres <- residuals(gnb)
ssdata$txeigres <- residuals(txlageig)
ssdata$txerres <- residuals(txerrorlag)

tm_shape(ssdata) + tm_polygons(col = "nbres", palette = "GnBu", title = "", border.col = "grey80", border.alpha = 0.5)

tm_shape(ssdata) + tm_polygons(col = "posres", palette = "GnBu", title = "", border.col = "grey80", border.alpha = 0.5)

tm_shape(ssdata) + tm_polygons(col = "txerres", palette = "GnBu", title = "", border.col = "grey80", border.alpha = 0.5)

tm_shape(ssdata) + tm_polygons(col = "txeigres", palette = "GnBu", title = "", border.col = "grey80", border.alpha = 0.5)

Trends

trend <- read_excel("C:/Users/xyaze/Desktop/trend.xlsx", sheet = "Sheet1")
trenddata <- merge(ssmap, trend)
trenddata <- st_as_sf(trenddata)
trendy <- read_excel("C:/Users/xyaze/Desktop/trend.xlsx", sheet = "Sheet4")
trendydata <- merge(ssmap, trendy)
trendydata <- st_as_sf(trendydata)
trend1 <- read_excel("C:/Users/xyaze/Desktop/trend.xlsx", sheet = "Sheet3")
trends <- merge(ssmap, trend1)
ggplot(trend1, aes(TX_CURR, fill = factor(State))) + geom_histogram() +  scale_fill_viridis_d(option = "viridis", direction = -1) + theme_bw() + facet_wrap(~YEAR) + labs(x = "TX_CURR", y = "Frequency") + theme(legend.position = "bottom", legend.title = element_blank(), legend.text = element_text(size = 0.5))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## Warning: Removed 103 rows containing non-finite values (`stat_bin()`).

trend1 <- read_excel("C:/Users/xyaze/Desktop/trend.xlsx", sheet = "Sheet3")
trends <- merge(ssmap, trend1)
tm_shape(trends) + tm_polygons('TX_CURR', title = "TX_CURR", palette = "-viridis", border.col = "white", border.alpha = 0.5) + tm_facets("YEAR", free.coords = FALSE) + tm_layout(panel.labels = c('2018', '2019', '2020', '2021', '2022'), legend.show = FALSE)

` ##